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a field of its own whithin applied mathematics, and is now essential to many parts of 
theoretical physics, from condensed matter to high energy. The fundamental results 
obtained so far rely mostly on the theory of random matrices in one dimension (the 
dimensionality of the spectrum, or equilibrium probability density). In the last few 
years, this theory has been extended to the case where the spectrum is two-dimensional, 
or even fractal, with dimensions between 1 and 2. In this article, we review these recent 
developments and indicate some physical problems where the theory can be applied. 
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1. Introduction 

During tlie second lialf of last century and continuing through the present, random 
matrix theory has grown from a special method of theoretical physics, meant to 
approximate energy levels of complex nuclei [H El [3l HJ El El [7] , into a vast mathematical 
theory with many different application in physics, computer and electrical engineering. 
Simply describing all the developments and methods currently employed in this context 
would result in a monography much more extensive than this review. Therefore, we will 
only briefly mention topics which are themselves very interesting, but lie beyond the 
scope of this work. 

The applications of random matrix theory (RMT) into physics have been extended 
from the original subject, spectra of heavy nuclei, to descriptions of large N SU{N) 
gauge theory [H [9], critical statistical models in two dimensions [lOl IHl [12] disordered 
electronic systems [131 HH [T51 [ISl 113 IB], quantum chromo-dynamics (QCD) [TSl [20], 
to name only a few. Non-physics applications range from communication theory [21] to 
stochastic processes out of equilibrium [221 [2S| and even more exotic topics [24J . 

A number of important results, both at theoretical and applied levels, were obtained 
from the connection between random matrices and orthogonal polynomials, especially 
in their weighted limit [23 12S1 13 123 [2H] ■ These works explored the relationship between 
the branch cuts of spectral (Riemann) curves of systems of differential equations and the 
support of limit measures for weighted orthogonal polynomials. Yet another interesting 
connection stemming from this approach is with the general (matrix) version of the 
Riemann-Hilbert problem with finite support [23 |30] . 

In [311 132] , it was showed that such relationships also hold for the class of normal 
random matrices. Unlike in previous works, for this ensemble, the support of the 
equilibrium distribution for the eigenvalues of matrices in the infinite-size limit, is two- 
dimensional, which allows to interpret it as a growing cluster in the plane. Thus, a direct 
relation to the class of models known as Laplacian Growth (both in the deterministic 
and stochastic formulations), was derived, with important consequences. In particular, 
this approach allowed to study formation of singularities in models of two-dimensional 
growth. Moreover, these results allowed to define a proper way of continuing the solution 
for singular Laplacian Growth, beyond the critical point. 

From the point of view of the dimensionality of the support for random matrix 
eigenvalues, it is possible to distinguish between 1-dimensional situations (which 
characterize 1 and 2-matrix models), and 2-dimensional situations, like in the case of 
normal random matrix theory. In fact, very recent results point to intermediate cases, 
where the support is a set of dimensional between 1 and 2. This situation is very similar 
to the description of disordered, interacting electrons in the plane, in the vicinity of 
the critical point which separates localized from de-localized behavior [15]. It is from 
the perspective of the dimensionality of support for equilibrium measure that we have 
organized this review. 

The paper is structured in the following way: after a brief summary of the main 
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concepts in Section [2| we explain the structure of normal random matrices in the limit 
of infinite size, in Section [3j This allows to connect with planar growth models, of 
which Laplacian (or harmonic) Growth is a main representative. The following two 
sections give a solid description of the physical (Section [i]) and mathematical (Section |5]) 
structure of harmonic growth. The discretized (or quantized) version of this problem is 
precisely given by normal random matrices, as we indicate in these sections. Next 
we present a general scheme for encoding shade functions in the plane into linear 
data, specifically into a linear bounded Hilbert space operator T with rank one self- 
commutator rank[T*,T] = 1. This line of research goes back to the perturbation and 
scattering theory of symmetric operators (M. G. Krein's phase shift function) and to 
studies related to singular integral operators with a Cauchy kernel type singularity. 
Multivariate refinements of the "quantization scheme" we outline in Chapter 5 lie at 
the foundations of both cyclic (co)homology of operator algebras and of free probability 
theory. In view of the scope and length of the present survey, we confine ourselves to 
only outline the surprising link between quadrature domains and such Hilbert space 
objects. 

We conclude with an application of the operator formalism to the description of 
boundary singular points that are characteristic to Laplacian growth evolution, and a 
brief overview of other related topics. 



2. Random Matrix Theory in ID 

2.1. The symmetry group ensembles and their physical realisations 

Following , we reproduce the standard introduction of the symmetry-groups ensemble 
of random matrices. The traditional ensembles (orthogonal, unitary and symplectic) 
were introduced mainly because of their significance with respect to symmetries of 
hamiltonian operators in physical theories: time-reversal and rotational invariance 
corresponds to the orthogonal ensemble (which, for Gaussian measures, is naturally 
abbreviated GOE), while time-reversal alone and rotational invariance alone correspond 
to the symplectic and unitary ensembles, respectively (GSE and GUE for Gaussian 
measures) . 

An invariant measure is defined for each of these ensembles, in the form 

dJl{M) = P{M)dfi{M) = Z-ie-T'-[^(^')l(i/x(M), (1) 

where M is a matrix from the ensemble, Z is a normalization factor (partition 
function), Tr[iy(M)] is invariant under the symmetries on the ensemble, and 
dii[M) is the appropriate flat measure for that ensemble: ni<j dMij for orthogonal, 
Y[i<j dRe Mij Yl-^j dim Mij for unitary, and Yl-^j dM^f HLi ^i<j dM^j ^ for symplectic 
(where each matrix element is an element of the real Klein group, Mjj = M^j'^ ■ 1 + 
J2k=i ^ij'^ ' '^k)- Correspondingly, to each of these ensembles, a parameter /3, indicating 
the number of independent real parameters necessary to describe the pair of values 
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Mij,Mji, is introduced, with values (3 = 1,2,4 for orthogonal, unitary and symplectic 
ensembles, respectively. 

The invariance under transformations from the appropriate symmetry group leads 
to the following simplification of the measure: for any of these ensembles, the generic 
matrix M can be diagonalized by a transformation M = U~^AU, with U from 
the same group, and A = diag(Ai, . . . , Atv)- The Jacobian of the transformation 
M — i> A,U (where U is said to carry the "angular" degrees of freedom of M) is 
J = YliKj ~ ■^jl^ — 1^(^)1'^) with A the Vandermonde determinant. The angular 
degrees of freedom can be integrated out (a trivial redefinition of the normalization 
factor), giving the simplified measure 

N N 

p(Ai, ...,X^)l[d\ = Z-ie^-^l^^^)! I A(A) f J] d\ (2) 



1=1 i=l 

r2 



For example, in the case of Gaussian measure W{M) = — M , the joint probability 
distribution function of eigenvalues, p, becames (up to normalization) 

TV 



p(Ai,...,A7v) =exp 



(3) 



5^A^ + /3^1og|A,-A, 

i=l i<j 

Clearly, this procedure is useful only if we are interested in computing expectation values 
of quantities which depend only of the distribution of eigenvalues, and not of the angular 
degrees of freedom. This is indeed the case for all situations of interest. 

The next standard transformation (which we discuss for the case of unitary 
ensemble, j3 = 2) that is performed on the measure uses the well-known property of 
Vandermonde determinant A (A) = det[A^~^]i<jj<Ar. Because of standard determinantal 
identities, this is equivalent with replacing each monomial A^~^ by a monic polynomial 
of the same order, Pj_i(Aj) = A^~^ + . . .. Finally, these polynomials may be chosen to 
be orthogonal with respect to the measure e^'^^\ giving for the p.d.f. of eigenvalues the 
expression 

p(Ai,...,A;v) = |det[P,_i(A,)e^(^»)/2]p, (4) 

which is simply the absolute value-squared of the wavefunction of the ground state 
for electrons in the external potential W . As we shall see, this kind o physical 
interpretation may be generalized to the case of matrix ensembles with two-dimensional 
support of eigenvalues. 



Generalizations of group ensembles Recently, various generalizations were proposed in 
order to extend the theory for ensembles of matrices which are not associated with 
symmetry groups. In particular, ensembles of matrices which may be reduced to a 
tridiagonal form (instead of standard diagonal) by a transformation which eliminates 
"angular" degrees of freedom, were introduced in [31]. As an interesting consequence, 
many results carry over to this case, while the parameter (3 is allowed to take any positive 
real value. 



2.2 Critical ensembles 
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2.2. Critical ensembles 

In this section we explain how, using properly chosen non-Gaussian measures, it is 
possible to construct ensembles of hermitian matrices (corresponding again to the 
unitary symmetry) which are in a sense, critical, i.e. for which a continuum limit 
(A^ oo) may be defined. The discussion relies on the formulation based on orthogonal 
polynomials indicated above, and it follows (at a more elementary level) the general 
theory of Saff and Totik [26]. 

2.2.1. General formalism Let dfj.{x) = e^^^^dx be a well-defined measure on the real 
axis, W{x) — > — oo as |x| oo, and p!i'\x) the corresponding family of orthogonal 
polynomials 

/>oo 

Pi'\x)Pj^\x)dflix) = 5nm. (5) 



Orthonormal functions are obtained through ipni^) = Pn{x)e^^^^^'^, which are 
orthogonal with respect to the fiat measure on M. We consider a deformation of this 
ensemble through a positive real parameter A > 1, so that dfi\{x) = e'^^^^^dx and 

Pi^\x)Pj^\x)dl^x{x) = 6nm. (6) 

Clearly, if W{x) is a monomial of degree A:, the deformation amounts to a simple rescaling 
p::{x) = X'/''p!^'\X'/'x). (7) 

The first non-trivial example is a quartic polynomial of the type 

W{x) = -(x^ + gx^), g>0, (8) 

for which the deformation in not a simple rescaling. In this case, it is possible to consider 
a special limit n — > oo, A oo, A nvc, where is a constant. As we will see, for a 
specific value of Vc, this limit yields a special asymptotic behavior of the orthonormal 
functions ipnix)- However, even for the simplest, trivial monomial (a Gaussian), which 
yields the Hermite polynomials, the asymptotic behavior of the orthogonal functions 
is non-trivial, in the sense that there are no known good approximations for the case 
re = 0(1). 

Generically, in this large n, A limit, we can ask where the wavefunction ipn{x) will 
reach its maximum value, in the saddle point approximation: 

max|a;|9^|?/'n(a;)| = 0, (9) 



giving 



so that 



log(x — ^i) — + W{x) 



i=l 



9 " 1 
- r^W'ix) = - y 



(10) 



:ir 



n ^ — ' X 

i=l 
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where ^j, z = 1, . . . , are the roots of the n^^ polynomial. 



Let 



1 " 1 

^{^) = ~ / 7 ) 



multiply [11] hj {^i — z) ^ and sum over i, and obtain 

W'iz)-W'{^i 



Lj\z)-r,W'{z)u{z) = --J2 



i=l 



(12) 



(13) 



Equation (13) can be solved in the large n limit by assuming that the roots will be 



distributed with density p(^) on some compact (possibly disconnected) set / C 
Defining 



R{z) = -- 



' c Jl 



W'{z) - W'{C) 



we obtain 



uj\z) - r,W\z)oj{z) + (l)'^(^) = 0. 



(14) 



(15) 



The proper solution of (13) (considering the behavior at oo of the function uj{z)), is 

r. 



uj{z) = ^ W\z) + ^{W'{z)Y-R{z) 



(16) 



and (since the function uj{z) is the Cauchy transform of the density p(x)), it gives us 
the asymptotic distribution of zeros as 

1 



p{x) 



[uj{x + iO) — uj{x — i( 



(17) 



Finally, to obtain the asymptotic form of wave functions ipnix), we can write 

n-HogM^)^ j PiO^ogix-Od^-'^ + Wix). (18) 

2.2.2. Continuum limit and integrable equations There are two related problems for 
the large n limit of deformed ensembles described in the previous section. The first 
is determination of the support of zeros /; the second is the scaling behavior of the 
orthogonal functions ipnix)- In general, the limiting support I may consist of several 
disconnected segments Ik, I = Ulzflk- In the simplest case, it is just one interval 
/ = [a, b] C M. In this section we indicate how to determine this support as well as the 
density p{x), and what this yields for the orthogonal functions. 



Let the function W{x) be a polynomial of even degree d. From (14) we see that 



R{z) is a polynomial of degree d — 2, and therefore solution (16) has generically 2{d — 1) 
branch points. Thus, the function uj{z) typically has d—1 branch cuts, which constitute 
the disconnected support of distribution p{z). 

We are interested in a special case, when li — 2 of these cuts degenerate into double 
points, and there is a single interval [a, b] which is the support of p{z). This special case 
is called critical and it provides new asymptotic limits for the orthogonal functions. We 
will also refer to this solution as the "single-cut" solution. 
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From the equation 

u{x + iO) - u{x - iO) = rcW'{x) 
we obtain for tlie single-cut solution 

rc^/{z -a){z - h) 



(19) 



uj{z) 



27r 



The large \z\ behavior of this function is known from the continuum limit of (12), and 
it implies the absence of regular terms in the Laurent expansion: 



uj{z) = -- + 0{z~^ 
so that we impose the conditions 







2n = —Tr 



a V(&-0(e-«) 



di. 



(20) 

(21) 
(22) 



Gaussian measure and the Hermite polynomials Let d = 2 and —W{x) = ax'^, a > 0. 
Then conditions ( |2Tp2| give a symmetric support [— 6, where 6^ = 2 /{arc). 



More generally, using the saddle point equation for |^/'„(x)| at x = a, 6 and (18), we 
conclude that 



^ JO 

Since the integrand behaves like rj'^-'^/'^^ we obtain 

nr, 



^nib + C) = V'n(&)exp 



2d - 1 



d-l/2 



(23) 



(24) 



We immediately conclude that for = 2 (Hermite polynomials), the asymptotic behavior 
is given by the Airy function, expz^/^. The full scaling is achieved by considering the 
region around the end-point b, of order ( = 0(?T,~^/'^^'^~^^).Then we obtain 



exp 



2d - 1 



(25) 



2.2.3. Scaled limits of orthogonal polynomials and equilibrium measures The 
distribution of eigenvalues investigated in the previous sections illustrates the general 
approach developed by Saff and Totik for holomorphic polynomials orthogonal on 
curves in the complex plane. We sketch here the more general result because of its 
relevance to the main topic of this review. 

Given a set S G C and a properly-defined measure on it w{z) = e~'^^^\ we construct 
the holomorphic orthogonal polynomials Pn{z), with respect to w. We then pose the 
question of finding the "extremal" measure (its support and density fi^), such that 
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the F-functional F{K) = logcap(fr) — J QduK, with cap{K) and ljk the capacity, 
respectively the equihbrium measure of the set K, is maximized by S^j. Furthermore, 
fiw satisfies energy and capacity constraints on Sw 

The remarkable fact noticed in is that if the extremal value F{Sy^) is 
approximated by the weighted monic polynomials Pn{z) as (| |w"P„| |^)^/" — ^ exp(— F^) 
(where we use the weak star norm), then the asymptotic zero distribution of P„ gives the 
support Syj. Hence, (17) may be interpreted as giving both the support of the extremal 
measure (labeled p in this formula), as well as its actual density. 

The extremal measure has the physical interpretation of the "smallest" equilibrium 
measure which gives a prescribed logarithmic potential at infinity. According to the 
concept of "sweeping" (or "balayage", see [2E]), the extremal measure is obtained as a 
limit of the process, under the constraints imposed on the total mass and energy of the 
measure. As we have shown in this chapter, for the case of ID measures, this extreme 
case is given by weighted limits of orthogonal polynomials. 



3. Random Matrix Theory in higher dimensions 

In this chapter, we show how to generalize the concepts of equilibrium measure, extremal 
measure, and their relations to orthogonal polynomials and ensembles of random 
matrices, in the case of two-dimensional support. The applications of this theory to 
planar growth processes will be discussed in the following two chapters. 



3.1. The Ginibre- Girko ensemble 

We begin with a brief discussion on the oldest and simplest ensemble of random matrices 
with planar support. The ensemble of complex, N x N random matrices with identical, 
independent, zero-mean Gaussian-distributed entries, was first studied by J. Ginibre in 
1965 [35j, and then it was generalized for non-zero mean Gaussian by Girko in 1985 [36] . 
Consider N x N random matrices with eigenvalues G C, and joint p.d.f. 

dPN ^ Y\ n Ai^(^fc)' (26) 

l<i<j<N l<k<N 

where fiNi^k) = e~'^^^''^^ dKeZkdlmZk. Then, in the large N limit, the measure 
jf Y2k ^(^ ~ ^k) converges weakly to the uniform measure on the unit disk. This 
is known as the Circular Law. If the exponent of the pure Gaussian is perturbed 
by a quadratic term, the result holds for a corresponding elliptical domain, giving 
the Elliptical Law. The same limiting curves (circular and elliptical) describe the 
graph of the distribution of real eigenvalues for Hermitian ensembles, with pure and 
perturbed Gaussian measures. In that case, the laws are known as Wigner-Dyson [H H] 
and Marchenko-Pastur, respectively (although the last one was originally derived for 
covariance matrices built from sparse regression matrices |37j). 

Extensions and exceptions from the circular and elliptical laws were found by 
relaxing the conditions of the theorems. In particular, deviations from uniformity for 
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angular statistics in the case of Gaussian measure were derived in [38], while the case 
of heavy-tails distributions was investigated in [2S1 SD] and subsequent publications. 

3.2. Normal matrix ensembles 

A special case of matrices with complex eigenvalues is given by normal matrices. A 
matrix M is called normal if it commutes with its Hermitian conjugate: [M, M'^] = 0, 
so that both M and M''' can be diagonalized simultaneously. The statistical weight of 
the normal matrix ensemble is given through a general potential W{M, M"^) |?T]: 

gitrW'(M,Mt)^^^^^_ (27) 

Here h is a parameter, and the measure of integration over normal matrices is 
induced by the fiat metric on the space of all complex matrices dcM, where dcM = 
Ylij dRe Mijdlm Mij. Using a standard procedure, one passes to the joint probability 
distribution of eigenvalues of normal matrices zi,...,zn, where N is the size of the 
matrix: 

1 ^ 

^lAivWPn^^"^^'^''^^^'^^- (28) 



Here d'^Zj = dxj dyj for zj = Xj + iyj, Ajv(z) = det{zj ^)i<i,j<N = Y[i>ji^i ~ is the 
Vandermonde determinant, and 



„ N 

rN = ^l \A^{z)\'lle'^^^^^'^^^d'z, (29) 

is a normalization factor, the partition function of the matrix model (a r-function). 
A particularly important special case arises if the potential W has the form 



W = -\z\' + V{z) + V{z), (30) 
where V{z) is a holomorphic function in a domain which includes the support of 



eigenvalues (see also a comment in the end of Section 3.4 about a proper definition 
of the ensemble with this potential). In this normal matrix ensemble gives the 

same distribution as a general complex matrix ensemble. A general complex matrix can 
be decomposed as M = U{Z + R)U\ where U and Z are unitary and diagonal matrices. 



respectively, and R is an upper triangular matrix. The distribution (28) holds for the 



elements of the diagonal matrix Z which are eigenvalues of M. Here we mostly focus 



on the special potential (30), and also assume that the field 

A{z) = d,V{z) (31) 
is a globally defined meromorphic function. 
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Figure 1. A support of eigenvalues consisting of four disconnected components (left). 
The distribution of eigenvalues for potential V{z) = — alog(l — z/f3) — jz. (right) 



3. 3. Droplets of eigenvalues 

In the large limit {h 0, Nh fixed), the eigenvalues of matrices from the ensemble 
densely occupy a connected domain D in the complex plane, or, in general, several 
disconnected domains. This set (called the support of eigenvalues) has sharp edges 
(Figure [T]). We refer to the connected components Da of the domain D as droplets. 

For algebraic domains (the definition follows) the eigenvalues are distributed with 
the density p = —^AW, where A = idzdg is the 2-D Laplace operator [3l]. For the 



potential (30) the density is uniform. The shape of the support of eigenvalues is the 
main subject of this chapter. For example, if the potential is Gaussian [55] . 

A{z) = 2t2Z, (32) 

the domain is an ellipse. If A has one simple pole, 

A{z) = -^-^ (33) 
z — p 

the droplet (under certain conditions discussed below) has the profile of an aircraft wing 
given by the Joukowsky map (Figure [T]). If A has one double pole (say, at infinity), 

A{z) = ?>t^z\ (34) 

the droplet is a hypotrochoid. If A has two or more simple poles, there may be more 
than one droplet. This support and density represent the equilibrium solution to an 
electrostatic problem, as we will indicate in a later section. 

3.4- Orthogonal polynomials and distribution of eigenvalues 
Define the exact A^-particle wave function (up to a phase), by 

^^(zi, ...,zm) = ^i=A^(z) e^f- hm-^.-^.). (35) 



3.5 Wavefunctions, recursions and integrable hierarchies 
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|2 



The joint probabihty distribution (28) is then equal to . . . iZ^q] 

Let the number of eigenvalues (particles) increase while the potential stays fixed. If 
the support of eigenvalues is simply-connected, its area grows as hN . One can describe 
the evolution of the domain through the density of particles 

Pn{z) = N j \^N{z,zi,Z2,...,ZN-i)\^d^zi...d^ZN-i, (36) 

where \E'7v is given by (35). 

We introduce a set of orthonormal one-particle functions on the complex plane as 
matrix elements of transitions between N and {N + l)-particle states: 

ip (z) f 

-^=== = y ^!N+l{z,Zi,Z2,...,ZN)'^N{zi,Z2,...,ZN)(f'Zi...(fzN (37) 

Then the rate of the density change is 

Pn+i{z) - Pn{z) = \tpN{z)\^- (38) 

The proof of this formula is based on the representation of the ipn through holomorphic 
biorthogonal polynomials Pn{z). Up to a phase 



^„(^) = e^^(^'^")p„(^), Pn{z) = U^z^ + ... (39) 

V 'Tn+l 

The polynomials Pn{z) are biorthogonal on the complex plane with the weight e^^^: 

j e'^'^Pn{z)Pj^)d^z = (40) 

The proof of these formulae is standard in the theory of orthogonal polynomials. 
Extension to the biorthogonal case adds no difficulties. 

We note that, with the choice of potential (30), the integral representation (40) has 
only a formal meaning, since the integral diverges unless the potential is Gaussian. A 
proper definition of the wave functions goes through recursive relations (53, 54) which 
follow from the integral representation. The same comment applies to the r-function 
(29). The wave function is not normalized everywhere in the complex plane. It may 
diverge at the poles of the vector potential field. 



3.5. Wavefunctions, recursions and integrable hierarchies 

In order to illustrate the mathematical connection between this theory and equivalent 
formulations which we present in Chapter 5, it is necessary to make a digression through 
the formalism of infinite, integrable hierarchies. In particular, we choose the case of the 
Kadomtsev-Petviashvilii (KP) hierarchy, and follow the notations in |42j . 

3.5.1. Pseudo-differential operators We denote by A the algebra constructed from 
differential polynomials of the type P = 9" -|-m„_2i9"~^ -|- . . . +uid + uo, where d = d/dz 
is a differential symbol with respect to some (complex) variable z, and uq, ui, . . . , u„_2 
(note: m„_i can be always set to zero) are generically smooth functions in z and (if 
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necessary) other variables ^1,^2, On this algebra, we define the ring of pseudo- 
differential operators TZ, consisting of (formal) operators defined by the infinite series 

i^=5^cfca^ d-'= / dz, (41) 
-00 

where coefficients are again smooth functions, and the negative powers in the expansions 
contain integral operators. For any such operator, we denote by L_|_ the purely 
differential part and by L_ the remainder of the series: 

n 

L+ = J2ckd', L = L+ + L_. (42) 



Let 



C = d + uod-^ + uid-^ + ... (43) 



be a pseudo-differential operator such that = d. Then, introducing the infinite set 
of times t = ti,t2, ■ . . , such that all coefficients Uk, Ck above are generically functions of 
t, the KP hierarchy has the form 

— = [4,>C], k^l,2,... (44) 

More explicitly, we note that the hierarchy consists of the differential equations satisfied 
by the coefficients of the operator £. As a consequence of the compatibility of all the 
equations in the hierarchy, we have the zero-curvature equations 

[d,^-C>l,dt^-Cl] = 0, Wtk,tp. (45) 

3.5.2. Level reductions The KP hierarchy contains many other known integrable 
hierarchies, particularly the KdV hierarchy, as reductions to a certain level n in the 
hierarchy. For example, assume that the operator C satisfies the constraint 

= 0, (46) 

i.e. it is the square root of a differential operator L of order 2: 

C = L^/^ L = + 2uo. (47) 

Then it follows that for all even powers n = 2m, L\ = so that [vC+,jC] = 0, so 
there is no dependence on the even times • • •• This sub-hierarchy is called level — 2 
KdV, because the first non-trivial zero-curvature equation of the hierarchy is the famous 
Korteweg-de Vries equation: 

3 dL 
£l = P^d^ + - [uod + duo] , ^ = [P, L] ^ ut, = 6uu, + u,,„ (48) 

where we have used Uq — u for clarity. 

This formulation of the KdV equation makes use of the notion of Lax pair L, P, 
which is central to the inverse scattering method for solving nonlinear integrable 
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differential equations The idea is quite physical: assume that the operators L, P act 
on a wavefunction t) such that 

= A^, ^ = P^, (49) 

where eigenvalues A form the spectrum of L. Then applying the Lax pair equation to 
the eigenvalue equation, we obtain dX/dt — 0, i.e. the evolution under these equations 
leaves the spectrum invariant. This allows to construct the initial state from the final 
state, hence the inverse scattering appellation. 

3.5.3. Tau functions and Baker-Akhiezer function At the level of systems of PDE, the 
r— function and the Baker- Akhiezer function are introduced, by analogy with the Lax 
par formulation indicated above, in the following way: 

Baker-Akhiezer function Consider the function il!{z, ti,t2, ■ ■ ■) satisfying 

Cip = zip, 1^ = Clip, Mk > 1. (50) 

Otk 

This is the Baker-Akhiezer function of the KP hierarchy. 

Fundamental property of the Baker-Akhiezer function 

Let = 1 + X]o° kid~^~^ be the "dressing" operator defined such that 
C = (f)d(f)~^ . Also, introduce the function g{z,ti,...) — exp[^^ ife^;*^]. Then 
the Baker-Akhiezer function satisfies: 

oo 

ip = k{z)g{z, ti, ...), k{z) = 1 + XI 



where k is the "scalar" analog of the dressing operator 0. 

Tau function Using the notation introduced above, we have the following property: 
There exists a function T{z,ti, . . .) such that 



r(^,ti,t2,---) r(^,ti,...) 
Now let us consider the generalized overlap function 



' ^) = '^N^ J "^N+liz, Zi, Z2,..., Zn)'^n+i{w, Zi, Z2, ■ ■ ■ , ZN)d'^Zi . . . d^ZN, 



iPn{z 



and expand for jz], l^l — > oo. We obtain 



ipN\z,w) — exp 



z'^wPdakp_ 

where Okp is the corresponding interior bi-harmonic moment. Therefore, we may regard 
the T-function and the scaled wavefunction introduced earlier as canonical objects 
describing an integrable hierarchy. This fact will be illustrated in more detail in the 
next section. 
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3. 6. Equations for the wave functions and the spectral curve 



In this section we specify the potential to be of the form (30). It is convenient to modify 
the exponential factor of the wave function. Namely, we define 

il^niz) = e-^+sn-)p„(^), and Xn{z) = es^(^)p„(z), (52) 

where the holomorphic functions Xn{z) are orthonormal in the complex plane with the 
weight e"'^'^/^. Like traditional orthogonal polynomials, the biorthogonal polynomials 
Pn (and the corresponding wave functions) obey a set of differential equations with 
respect to the argument z, and recurrence relations with respect to the degree n. Similar 
equations for two- matrix models are discussed in numerous papers (see, e.g., [13]). 

We introduce the L-operator (the Lax operator) as multiplication by z in the basis 

X.n- 

LnmXm{z) = ZXn{z) (53) 

(summation over repeated indices is implied). Obviously, L is a lower triangular matrix 
with one adjacent upper diagonal, Lnm = Oasm>n + l. Similarly, the differentiation 
dz is represented by an upper triangular matrix with one adjacent lower diagonal. 
Integrating by parts the matrix elements of the 5^, one finds: 

{L'^)nmX'm = ^T'dzXn, (54) 

where is the Hermitian conjugate operator. 

The matrix elements of L'l" are (!/''')„„ = Lmn = ^(-^nm) + / ^t^^ Pm{z)dzPniz)(P z, 
where the last term is a lower triangular matrix. The latter can be written through 
negative powers of the Lax operator. Writing dz\ogPn{z) = f + J2k>i'^k{n')z~'', one 
represents in the form 

= A{L) + {hn)L-^ + (55) 

fc>i 

where w'-'^-* and [hn] are diagonal matrices with elements fi^^ and {hn). The coefficients 
v'^'^ are determined by the condition that lower triangular matrix elements of A{Lnm) 
are cancelled. 

In order to emphasize the structure of the operator L, we write it in the basis of 
the shift operator such that wfn = fn+iw for any sequence /„. Acting on the wave 
function, we have: 

WXn = Xn+l- 

In the ^-representation, the operators L, L'^ acquire the form 

L = r„w + ^MW^-fc^ = w-^rn + Y^w^u^^\ (56) 

k>0 k>0 

X The shift operator w has no inverse. Below is understood as a shift to the left defined as 
w^^w = 1. Same is applied to the operator L^^ . To avoid a possible confusion, we emphasize that 
although Xn is a right-hand eigenvector of L, it is not a right-hand eigenvector of . 
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Clearly, acting on Xni we have the commutation relation ("the string equation") 

[L, Lt] = h. (57) 



This is the compatibility condition of Eqs. (53) and (54). 



Equations (56) and (57) completely determine the coefficients v, 



(k) 

n , 



and u 



(k) 



The ffist one connects the coefficients to the parameters of the potential. The second 
equation is used to determine how the coefficients vii 
particular, the diagonal part of it reads 



(k) 



Tn and Un'^ evolve with n. In 



nU 



EE 

fc>i p=i 



u, 



(k) 1 2 
n+p I • 



(5J 



Moreover, we note that all the coefficients can be expressed through the r- function (29) 



and its derivatives with respect to parameters of the potential. This representation is 
particularly simple for r„: = r„r^_^]^r„+2. 



3.6.1. Finite dimensional reductions If the vector potential A{z) is a rational function, 
the coefficients Un"^ are not all independent. The number of independent coefficients 
equals the number of independent parameters of the potential. For example, if the 



holomorphic part of the potential, V{z)., is a polynomial of degree d, the series (56) are 
truncated aX k = d — 1. 



In this case the semi- infinite system of linear equations (54) and the recurrence 



relations (53) can be cast in the form of a set of finite dimensional equations whose 



coefficients are rational functions of 2;, one system for every n > 0. The system 
of differential equations generalizes the Cristoffel-Daurboux second order differential 
equation valid for orthogonal polynomials. This fact has been observed in recent papers 
[m HH] for biorthogonal polynomials emerging in the Hermitian two-matrix model 
with a polynomial potential. It is applicable to our case (holomorphic biorthogonal 
polynomials) as well. 

In a more general case, when A{z) is a general rational function with d — 1 poles 



(counting multiplicities), the series (56) is not truncated. However, L can be represented 
as a "ratio". 



where the operators K12, Mi 2 are polynomials in w: 



d-2 



j=Q 
d-2 



Ml 



w 



-'+Ec^;. 

j=0 



Mo 



d-l 
j=0 



W-' 



d-l 

+ E^ 

j=0 



These operators obey the relation 
K1M2 = K2M1. 



(59) 



(60) 



(61) 



(62) 
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It can be proven that the pair of operators Mi 2 is uniquely determined by K12 and vice 



versa. We note that the reduction (59) is a difference analog of the "rational" reductions 



of the Kadomtsev-Petviashvili integrable hierarchy considered in [46] . 



The linear problems ( 53 ) , ([54j) acquire the form 

{K2X)n = Z {K,x)n , {Mlx)n = hd,{Mlx)n. (63) 

These equations are of finite order (namely, of order d), i.e., they connect values of Xn 
on d + 1 subsequent sites of the lattice. 

The semi-infinite set {xo, Xii • • •} is then a "bundle" of d-dimensional vectors 

X{n) = (Xn, Xn+1, Xn+d~l] 



(the index t means transposition, so x is a column vector). The dimension of the vector 
is the number of poles of A{z) plus one. Each vector obeys a closed ci-dimensional linear 
differential equation 

hd,x{n) = Cn{z)x{n), (64) 

where the dxd matrix Cn is a "projection" of the operator L"!" onto the n-th rf-dimensional 
space. Matrix elements of the are rational functions of z having the same poles as 
A{z) and also a pole at the point A{oo). (If A{z) is a polynomial, all these poles 
accumulate to a multiple pole at infinity). 

We briefly describe the procedure of constructing the finite dimensional matrix 



differential equation. We use the first linear problem in (63) to represent the shift 
operator as a d x d matrix yVn{z) with 2;-dependent coefficients: 

Wn{z)x{n) = x{n+l). (65) 
This is nothing else than rewriting the scalar linear problem in the matrix form. Then 



the matrix Wniz) is to be substituted into the second equation of (63 ) to determine C 



(examples follow). The entries of yVn{z) and Cn{z) obey the Schlesinger equation, which 



follows from compatibility of (64) and mb 



hd^Wn = Cn+lWn - WnCn. (66) 

This procedure has been realized explicitly for polynomial potentials in recent 
papers [111113]. We will work it out in detail for our three examples: x("-) = iXn, Xn+iY 



hypotrochoid (34). 



for the ellipse (32) and the aircraft wing (33) and xi^^) = iXn, Xn+i,Xn+2)* for the 



3.7. Spectral curve 

According to the general theory of linear differential equations, the semiclassical (WKB) 



asymptotics of solutions to Eq. (64), as h — > 0, is found by solving the eigenvalue 
problem for the matrix jCn{z) |17j. More precisely, the basic object of the WKB approach 
is the spectral curve [17] of the matrix Cn, which is defined, for every integer n > 0, by 
the secular equation det(£„(2;) — z) = (here z means z ■ 1, where 1 is the unit d x d 
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matrix). It is clear that the left hand side of the secular equation is a polynomial in z 
of degree d. We define the spectral curve by an equivalent equation 

fn{z,z) = a{z)&Qi{lin{z) - z) (67) 

where the factor 0(2) is added to make fn{z,z) a polynomial in z as well. The factor 
a{z) then has zeros at the points where poles of the matrix function C{z) are located. It 
does not depend on n. We will soon see that the degree of the polynomial a{z) is equal 
to d. Assume that all poles of A{z) are simple, then zeros of the a{z) are just the d — 1 
poles of A{z) and another simple zero at the point A{oo). Therefore, we conclude that 
the matrix Cn{z) is rather special. For a general d x d matrix function with the same d 
poles, the factor a{z) would be of degree d'^. 

Note that the matrix Cn{z) — z enters the differential equation 

nd,\±[n)\^ = l{n){C^{z) - z)l{n) (68) 
for the squared amplitude |'?/'(?T,)p = ilj'{n)ilj{n) = e^^|x(^)P of the vectors ilj{n) built 



from the orthonormal wave functions (39) 



The equation of the curve can be interpreted as a "resultant" of the non- 
commutative polynomials K2 — zKi and M\ — zMI (cf. [44J). Indeed, the point {z^z) 
belongs to the curve if and only if the linear system 



{K2c)k = z{Kic)k n — d < k < n 



(69) 



{Mlc)k = z{Mlc)k n<k <n + d-l 

has non-trivial solutions. The system contains 2d equations for 2d variables 
Cn-d , ■ ■ ■ ,Cn+d-i- Vanishing of the 2d x 2d determinant yields the equation of the 
spectral curve. Below we use this method to find the equation of the curve in the 
examples. It appears to be much easier than the determination of the matrix £n{z)- 



The spectral curve (67) possesses an important property: it admits an 
antiholomorphic involution. In the coordinates z,z the involution reads {z,z) 1— > {z,z). 
This simply means that the secular equation det(£„(5) — z) = for the matrix 
jCn{z) = Cn{z) defines the same curve. Therefore, the polynomial takes real values 
for z = z: 



Uz,z) = U{z,z). (70) 
Points of the real section of the curve {z = z) are fixed points of the involution. 



The curve (67) was discussed in recent papers jlH SS] in the context of Hermitian 
two-matrix models with polynomial potentials. The dual realizations of the curve 
pointed out in correspond to the antiholomorphic involution in our case. The 
involution can be proven along the lines of these works. The proof is rather technical 
and we omit it, restricting ourselves to the examples below. We simply note that the 
involution relies on the fact that the squared modulus of the wave function is real. 

We will give a concrete example for the construction of the spectral curve, after a 
brief but necessary detour through the continuum limit of this problem. 
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3.7.1. Schwarz function The polynomial fn{z,z) can be factorized in two ways: 

Uz, z) = a{z){z - S^:\z)) ...{z- Si'\z)), (71) 

where Sn\z) are eigenvalues of the matrix Cn{z), or 

Uz,z) = - SW(-z)) ...{z- Si'\-z)), (72) 

where Sn ^ (z) are eigenvalues of the matrix £„(^). One may understand them as different 
branches of a multivalued function S{z) (respectively, S{z)) on the plane (here we do 
not indicate the dependence on n, for simplicity of the notation). It then follows that 
S{z) and S{z) are mutually inverse functions: 

S{S{z)) = z. (73) 

An algebraic function with this property is called the Schwarz function. By the 
equation f{z, S{z)) = 0, it defines a complex curve with an antiholomorphic involution. 
An upper bound for genus of this curve is g = (rf — 1)^, where d is the number of branches 
of the Schwarz function. The real section of this curve is a set of all fixed points of the 
involution. It consists of a number of contours on the plane (and possibly a number of 
isolated points, if the curve is not smooth). The structure of this set is known to be 
complicated. Depending on coefficients of the polynomial, the number of disconnected 
contours in the real section may vary from to (? + 1. If the contours divide the complex 
curve into two disconnected "halves" , or sides (related by the involution) , then the curve 
can be realized as the Schottky double [48j of one of these sides. Each side is a Riemann 
surface with a boundary. 



Let us come back to equation (64). It has d independent solutions. They are 
functions on the spectral curve. One of them is a physical solution corresponding to 
biorthogonal polynomials. The physical solution defines the "physical sheet" of the 
curve. 

The Schwarz function on the physical sheet is a particular root, say Sn\z), of the 



polynomial fn{z,z) (see (71)). It follows from (55) that this root is selected by the 



requirement that it has the same poles and residues as the potential A. 

3.7.2. The Schottky double The Schwarz function describes more than just the 
boundary of clusters of eigenvalues. Together with other sheets it defines a Riemann 
surface. If the potential A{z) is meromorphic, the Schwarz function is an algebraic 
function. It satisfies a polynomial equation f{z, S{z)) = 0. 

The function f{z,z), where z and z are treated as two independent complex 



arguments, defines a Riemann surface with antiholomorphic involution (70). If the 
involution divides the surface into two disconnected parts, as explained above, the 
Riemann surface is the Schottky double [M] of one of these parts. 



There are two complementary ways to describe this surface. One is through the 



algebraic covering (71 72 ) . Among d sheets we distinguish a physical sheet. The physical 
sheet is selected by the condition that the differential S{z)dz has the same poles and 
residues as the differential of the potential A{z)dz. It may happen that the condition 
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Figure 2. The Schottky double. A Riemann surface with boundaries along the 
droplets (a front side) is glued to its mirror image (a back side). 

z = S^'^\z) defines a planar curve (or several curves, or a set of isolated points) for 
branches other than the physical one. We refer to the interior of these planar curves as 
virtual (or unphysical) droplets situated on sheets other than physical. 

Another way emphasizes the antiholomorphic involution. Consider a meromorphic 
function h{z) defined on a Riemann surface with boundaries. We call this surface the 
front side. The Schwarz reflection principle extends any meromorphic function on the 
front side to a meromorphic function on the Riemann surface without boundaries. This 
is done by adding another copy of the Riemann surface with boundaries (a back side), 
glued to the front side along the boundaries. Figure |2} The value of the function h on 
the mirror point on the back side is h{S{z)). The copies are glued along the boundaries: 
h{z) = h{S{z)) if the point z belongs to the boundary. The same extension rule apphes 
to differentials. Having a meromorphic differential h{z)dz on the front side, one extends 
it to a meromorphic differential h{S{z))dS{z) on the back side. 

This definition can be applied to the Schwarz function itself. We say that the 
Schwarz function on the double is S{z) if the point is on the front side, and z if the 
point belongs to the back side (here we understand S{z) as a function defined on the 
complex curve, not just on the physical sheet). 

The number of sheets of the curve is the number of poles (counted with their 
multiplicity) of the function A{z) plus one. Indeed, poles of A are poles of the Schwarz 
function on the front side of the double. On the back side, there is also a pole at infinity. 
Since S{z = oo) = A{oo), we have S{z = A{oo)) = oo. Therefore, the factor a{z) is a 
polynomial with zeros at the poles of A{z) and at A{oo), and 

d = number of sheets = number of poles of A + 1. 

The front and back sides meet at planar curves z = S{z). These curves are boundaries 
of the droplets. We repeat that not all droplets are physical. Some of them may belong 
to unphysical sheets. Figure [3] 
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Figure 3. Physical and unphysical droplets on a torus. The physical sheet (shaded) 
meets the unphysical sheet along the cuts. The cut situated inside the unphysical 
droplet appears on the physical sheet. The boundaries of the droplets (physical and 
virtual) belong to different sheets. This torus is the Riemann surface corresponding to 
the ensemble with the potential V{z) = — alog(l — z//3) — jz. 

Boundaries of droplets, physical and virtual, form a subset of the a-cycles on the 
curve. Their number cannot exceed the genus of the curve plus one: 

number of droplets < g + I. 

The sheets meet along cuts located inside droplets. The cuts that belong to physical 
droplets show up on unphysical sheets. On the other hand, some cuts show up on the 
physical sheet (Figure |3|. They correspond to droplets situated on unphysical sheets. 

The Riemann-Hurwitz theorem computes the genus of the curve as 

g = half the number of branching points — d + 1. 

With the help of the Stokes formula, the numbers are identified with areas 
of the droplets: [t'al = 2th Id ^"^^^ ^ nondegenerate curve, these numbers are not 
necessarily positive. Negative numbers correspond to droplets located on unphysical 
sheets. In this case, {ua} do not correspond to the number of eigenvalues located inside 
each droplet, as it is the case for algebraic domains, when all filling numbers are positive. 

3. 7. 3. Degeneration of the spectral curve Degeneration of the complex curve gives the 
most interesting physical aspects of growth. There are several levels of degeneration. 
We briefly discuss them below. 

Algebraic domains and double points A special case occurs when the Schwarz function 
on the physical sheet is meromorphic. It has no other singularities than poles of A. 
This is the case of algebraic domains . They appear in the semiclassical case. This 
situation occurs if cuts on the physical sheet, situated outside physical droplets, shrink 
to points, i.e., two or more branching points merge. Then the physical sheet meets other 
sheets along cuts situated inside physical droplets only and also at some points on their 
exterior [double points). In this case the Riemann surface degenerates. The genus is 
given by the number of physical droplets only. The filling factors are all positive. 



3.7 Spectral curve 



23 




Figure 4. Degenerate torus corresponds to the algebraic domain for the Joukowsky 
map. 



In the case of algebraic domains, the physical branch of the Schwarz function 
is a well-defined meromorphic function. Analytic continuations of z from different 
disconnected parts of the boundary give the same result. In this case, the Schwarz 
function can be written through the Cauchy transform of the physical droplets: 



Although algebraic domains occur in physical problems such as Laplacian growth, 
their semiclassical evolution is limited. Almost all algebraic domains will be broken in 
a growth process. Within a finite time (the area of the domain) they degenerate further 
into critical curves. The Gaussian potential (the Ginibre-Girko ensemble), which leads 
to a single droplet of the form of an ellipse is a known exception. 

Critical degenerate curves Algebraic domains appear as a result of merging of simple 
branching points on the physical sheet. The double points are located outside physical 
droplets. Remaining branching points belong to the interior of physical droplets. 
Initially, they survive in the degeneration process. However, as known in the theory 
of Laplacian growth, the process necessarily leads to a further degeneration. Sooner or 
later, at least one of the interior branching points merges with one of the double points 
in the exterior. Curves degenerated in this manner are called critical. For the genus one 
and three this degeneration is discussed below. 

Since interior branching points can only merge with exterior branching points on 
the boundary of the droplet, the boundary develops a cusp, characterized by a pair p, q 
of mutually prime integers. In local coordinates around such a cusp, the curve looks like 
y'i. The fact that the growth of algebraic domains always leads to critical curves 
is known in the theory of Laplacian growth as finite time singularities. 

The degeneration process seems to be a feature of the semiclassical approximation. 
Curves treated beyond this approximation never degenerate. 

3.7.4. Example: genus one curve The potential is V{z) = — alog(l — z/(3) — 
7^, A{z) = ~ 7- There is one pole at 2; = /? on the first (physical) sheet. 

At 2; = cxD on the first sheet S{z) — > —7 + Therefore, the Schwarz function has 

another pole at the point —7 on another sheet. All the poles are simple. According to 




(74) 
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the general arguments of Sec. 3.7.2 the number of sheets is 2, the number of branching 



points is 4. The genus is 1. The curve has the form 

f{z, z) = z^T?' + k\z^z + kxzz"^ + k2Z^ + k2Z^ + k^^zz + k^z + k^z + /i = 0. 

The points at infinity and —7 belong to the second sheet of the algebraic covering. 
Summing up, 



(-7 + ^) as z^ooi, 
as Z--72, 
(/3 - f ) as z ^ 002. 



where, by 1 and 2 we indicate the sheets. 

Poles and residues of the Schwarz function determine all the coefficients of the curve 



z) = a{z){z - S^^\z)){z - 5(2) (z)) = a{z){z - S^^\z)){z - S^^\z)) except one. The 
behavior at 00 of z,z gives fci = 7 — /5, = —iP- Hereafter we choose the origin by 
setting 7 = 0. The equation of the curve then reads fn{z, z) = 0, where /„(z, z) is given 
by 

z'^z'^—z'^z/3—zz'^P+(\P\'^ + a + a — nh) zz+z/3{nh—a)+zP{nh—a)+hn(75) 

The free term hn is to be determined by filling factors of the two droplets z/i and 
i>2 = n — ui. A detailed analysis shows that the droplets belong to different sheets 
(Figure [3]). Therefore, z/2 is negative. 

A boundary of a physical droplet is given by the equation z = S^^'^{z) (Figure 
The second droplet belongs to the unphysical sheet. Its boundary is given hy z = S^'^\z). 
The explicit form of both branches is 

c(i 2) _ 1 75 _ /3(^^ -a) + {a + a- nh)z ^ ^/ {z - zi){z - Z2){z - z:i){z - Zj) 
r 2{z-(3)z ' 

where the branching points Zi depend on hn- 

If the filling factor of the physical droplet is equal to ra, the cut inside the unphysical 
droplet is of the order of \/h. Although it never vanishes, it shrinks to a double point 
^3 = 2:4 = 2;* in a semiclassical limit. The sheets meet at the double point z^ rather than 
along the cut: {z — Zi){z — Z2) {z — 2:3) {z — Z4) {z — z^^{z — Z\){z — ^2)- In this 
case, genus of the curve reduces to zero and the exterior of the physical droplet becomes 
an algebraic domain. This condition determines Zi, and also the position of the double 
point (Figure [4]). The double point is a saddle point for the level curves of f{z,z). If 
all the parameters are real, the double point is stable in x-direction and unstable in 
y-direction. 

If this solution is chosen, the exterior of the physical droplet can be mapped to the 
exterior of the unit disk by the Joukowsky map 

u 

z(w) = rw + uo -\ , \w\ > 1, \a\ < 1. (76) 

w — a 
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The inverse map is given by the branch Wi{z) (such that wi 
double valued function 



oo as 2; ^ 00) of the 



1 

2r 



z - Uo + ar ± >y{z - zi){z - Z2) , -21,2 = Uq + ar ^ 2\/r{u + au^. 



The function 



z\w 



rw ^ + Mo + 



u 



1 - (77) 

is a meromorphic function of w with two simple poles at w = and w = d'^. Treated 
as a function of z, it covers the 2;-plane twice. Two branches of the Schwarz function 
are S^^''^\z) = z{wil{z)). On the physical sheet, S^^\z) = z{wi{z)) is the analytic 
continuation of z away from the boundary. This function is meromorphic outside the 
droplet. Apart from a cut between the branching points 2:1^2; the sheets also meet at 
the double point z^: = —7 + a~^re^^'^, where S^^\z^:) = S'''^\z^:), cp = arg(ar + jq^)- 

Analyzing singularities of the Schwarz function, one connects parameters of the 
conformal map with the deformation parameters: 
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A critical degeneration occurs when the double point merges with a branching point 
located inside the droplet (z* = Z2) to form a triple point z^^. This may happen on the 
boundary only. At this point, the boundary has a (2, 3) cusp. In local coordinates, it 
is ~ j/^. This is a critical point of the conformal map: w'{z^:^) = 00. A critical point 
inevitably results from the evolution at some finite critical area. 

A direct way to obtain the complex curve from the conformal map is the following. 



First, rewrite (76) and (77) as 



z — uq + ar 
z — Uq + dr 



rw + a{z + 'y)w ^ 
rw^^ + d{z + 7)iy, 



(79) 



and treat w and 1/w as independent variables. Then impose the condition w -w ^ 
One obtains 



1. 



det 



z — uq + ar a{z + 7) 
z ~ Uq + dr r 



det 



r a{z + 7) 
d{z + 7) r 



This gives the equation of the curve and in particular h, in terms of u, Uq, r, a and 
eventually through the deformation parameters a, P, 7 and t. 

The semiclassical analysis gives a guidance for the form of the recurrence relations. 



Let us use an ansatz for the L-operator, which resembles the conformal map (|76|): 

L 



VnW + + (W - an) 
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so that 

where w is the shift operator n — n + 1. 

Now we follow the procedure of the previous section. Since the potential has only 



(80) 
(81) 



one pole, Cn can be cast into 2x2 matrix form. Let us apply the lines (80, 81) to an 



[Z + Tn-lttn-l -U'n )Cn 



2) 



eigenvector (c„, c„+i) of a yet unknown operator and set the eigenvalue to be z: 

(0)n _ I f , - \ 

[Z + Tnan - U''n+l)^ri = a„+i (5 + 7„+i)Cn+l 

We have defined 7n = ^ — The equations are compatible if c„_i and c„+i found 
through Cn differ by the shift n ^ n + 2. We have 

,(0) 



Cn+l 



Cn-l 



d. 



^det 



z + rn-ittn-i-uli' a„-i(2; + 7n-i) 



det 



a„+i(5 + 7„+i) z + rndn-Unli 



(0) 



On 
dr,. 



where 



dr, = det 



This yields the curve 



V„ ■ V, 



n+l 



an+i{z + 7„+lj 



dndn+1- 



(83) 
(84) 

(85) 
(86) 



Comparing the two forms of the curve (75) and (86), we obtain the conservation laws 
of growth: 



1 = ln 



Ur. 



- u: 



^7) 



p 



a. 



(0) 

(0) I U^^ittnCln+l 



n+l 



nh-a = r„r„+i 



1 ~ O-ndn+l 



^9) 



They are the quantum version of (78). 



3.8. Continuum limit and conformal maps 



The geometrical meaning of the complex curve (67) is straightforward: at fixed shape 



parameters tk and area parameter h, increasing n yields growing domains that represent 
the support of the corresponding n x n model. A remarkable feature of this process is 
that it preserves the external harmonic moments of the domain D„, 

1 r d^z 



tk{n) =tk{n-l), tk{n) 



k> 1. 



(90) 
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The only harmonic moment which changes in this process is the normahzed area 
to = ^ J d'^z, and it increases in increments of h (hence the meaning of h as quantum 
of area). We may say that the growth of the NRM ensemble consists of increasing the 
area of the domain by multiples of h, while preserving all the other external harmonic 
moments. The continuum version of this process, known as Laplacian Growth, is a 
famous problem of complex analysis. It arises in the two-dimensional hydrodynamics 
of two non-mixing fluids, one inviscid and the other viscous, upon neglecting the effects 
of surface tension, where it is known as the Hele-Shaw problem. The following chapters 
discuss this classical problem in great detail. 

As we will see, Laplacian growth can be restated simply as a problem of finding 
the uniform equilibrium measure, subject to constraints on the total mass, and the 
asymptotic expansion of the logarithmic potential at infinity. As long as a classical 
solution exists, the machinery of NRM does not seem necessary. However, Laplacian 
growth (as a class of processes), is characterized by finite-time singularities. In that case, 
the only way to reformulate the problem is similar to the Saff-Totik approach to the 
extremal measure, and is deeply related to weighted limits of orthogonal polynomials in 
the complex plane. 

4. Laplacian Growth 

4.I. Introduction 

Laplacian growth (LG) is defined as the motion of a planar domain, whose boundary 
velocity is a gradient of the Green function of the same domain (also called a harmonic 
measure). This deceivingly simple process appears to be connected to an impressive 
number of non-trivial physical and mathematical problems |l9l[50]. As a highly unstable, 
dissipative, non-equilibrium, and nonlinear phenomenon, it is famous for producing 
different universal patterns [HH [52] . 

Numerous non-equilibrium physical processes of apparently different nature are 
examples of Laplacian growth: viscous fingering |[51J, slow freezing of fluids (Stefan 
problem) [53], growth of snowfiakes |5l], crystal growth, amorphous solidification [55] , 
electrodeposition [HB], bacterial colony growth [57J, diffusion-limited aggregation (DLA) 
[58] . motion of a charged surface in liquid Helium [59], and secondary petroleum 
production [60], to name just a few. 

A major consequence from the current development of the subject is a discovery of 
a new and unexpectedly fruitful mathematical structure, which is capable to predict and 
explain key physical observations in regimes, totally inaccessible by any other available 
mathematical method. 

The first section of this chapter is a brief history of physics covered by the Laplacian 
growth. The second section addresses in detail the exact time- dependent solutions of 
the Laplacian Growth Equation, and the last is a detailed presentation of the analytic 
and algebraic-geometric structure of Laplacian growth. 
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4-2. Physical background 

Darcy's law In 1856, while completing a hydrological study for the city of Dijon, H. 
Darcy noticed that the rate of flow (volume per unit time) Q through a given cross- 
section, to; is (a) proportional to Iq, (b) inversely proportional to the length, L, taken 
between positions of efflux and influx, and (c) linearly proportional to pressure difference, 
Ap, taken between the same two levels. In short, 

Q = ~^P, (91) 

where is a positive constant. As one can see, Darcy's observation coincides with Ohm's 
law, upon identifying Q, Ap, and k as the total current through the cross-section to, 
the electric potential difference, and the electrical conductivity, respectively. Rewriting 



(91) in a differential form, as for Ohm's law, we obtain 

V = -kVp, (92) 

where v is the velocity vector fleld of fluid particles, properly coarse-grained to assure 
its smoothness over inflnitesimally small volumes. Here the kinetic coefficient, k, (the 



same as in (92)), is called a (hydraulic) conductivity and can depend on position. The 
equation (92) constitutes the Darcy's law in a differential form. For homogenous k, 

V = Wi-kp), (93) 

Darcy's law merely states that a flow through uniform porous media (sand in Darcy's 
experiments) is purely potential (no vortices), where the pressure fleld, p, is a velocity 
potential up to a constant factor. Assuming constant k and the fluid incompressible, 
V ■ V = 0, we flnd that pressure p is a harmonic function, 

= (94) 

As seen from purely dimensional considerations, the conductivity k equals 

k = C-, (95) 

where d is the average linear size of a pore in cross-section, /i is the dynamical viscosity of 
the fluid under consideration, and the dimensionless coefficient, C, is usually small and 
media-dependent. (It is of the order of the density of voids in a given porous medium). 



It follows from ( |92| and (95), that if n is negligibly small (an almost inviscid liquid), 
pressure gradients are also negligibly small, regardless of how fast fluid moves (but still 
much slower than the velocity of sound in this liquid in order to assure incompressibility 
assumed earlier). 



Laplacian growth in porous media Assume that a fluid with a viscosity fii occupying 
a domain -Di(t) at the moment t pushes another fluid with a viscosity /i2 occupying the 
domain D2{t) at the same time t through a uniform porous media. Then the Laplace 
equation will hold for both pressures pi and p2 corresponding to domains Di and D2 
respectively: 

V^Pi = in Di{t), (96) 
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where i = 1,2. At the interface T{t), where two fluids meet (but do not mix), their 
normal velocities coincide because of continuity and equal to the normal component Vn 
of the velocity of the boundary, T(t) = dDi = —dD2: 

viU = V2U = K atr(t). (97) 

The pressure field p at the interface T{t) (by the Laplace law) has a jump equal to the 
mean local curvature k multiplied by the surface tension a: 

Pi—P2 = o'i^ at r(t). (98) 

Unless the local curvature is very high, this surface tension correction is usually very 
small, and so is often neglected. If to supplement the last three equations by boundary 
conditions at external walls or /and at infinity (they may include sources/sinks of fiuids 
either extended or point-like), then the free boundary problem of finding r(t) by initially 
given Di and D2 is completely formulated. 



The process described by (96 97, 98) is typical for various geophysical systems, for 
instance for petroleum production, where a less viscous fiuid (usually water) pushes a 
much more viscous one (oil) toward production wells. This process is very unstable and 
most initially smooth water/oil fronts will quickly break down and become fragmented. 

The Hele-Shaw cell In 1898, H.S. Hele-Shaw proposed an interesting way to observe 
and study two-dimensional fiuid fiows by using two closely-placed parallel glass plates 
with a gap between them occupied by the fiuid under consideration [61]. This simple 
device appears to be very useful in various investigations and is now called a Hele-Shaw 
cell after its inventor. Remarkably, a viscous fiuid, governed in 3D by the Stokes law, 

= Vp, (99) 

after being trapped in a gap of a width b, between the plates of a Hele-Shaw cell, obeys 



Darcy's law (92) with a conductivity equal to k = b'^/{12fj,). The derivation of the 
formula 

v = --^Vp, (100) 
which is to be understood as a 2D vector field in a plane parallel to the Hele-Shaw 



cell plates, is rather trivial and results from the averaging of (99) over the dimension 
perpendicular to the plates [621 [63]). Thus, displacement of viscous fiuid by the (almost) 
inviscid one in a Hele-Shaw cell became a major experimental tool to investigate a 
2D Laplacian growth. Various versions of 2D Laplacian growth in a Hele-Shaw cell 
corresponding to different geometries are shown in Figure |5j 



Idealized Laplacian growth In 1945, Polubarinova-Kochina [64] and Galin [65] 
simultaneously, but independently, derived a nonlinear integro-differential equation for 
an oil/water interface in 2D Laplacian growth, after neglecting surface tension, a, and 
water viscosity, fiwater = 0. Assuming for simplicity a singly connected oil bubble. 
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occupying a domain D{t) surrounded by water and having a sink at the origin, G D{t), 
we wiU obtain this equation starting from the system 

V^p = p in D{t), 

p = at the interface, T{t) = dD{t) (101) 
Vn = —dnP at the interface, r(t), 

where p and 9„ are density of sources and the normal derivative respectively. This 



system is a reduction of (96, 97) after simplifications mentioned above and using the 



fact that the normal boundary velocity, Vn, equals to the normal components of the 



fluid velocity at the boundary, which is —dnP by virtue of the Darcy law (92). Here and 
below the conductivity k is scaled to one. The density of sources, p, in this case equals 
p{z) = —5'^{z), which corresponds to a sink of unit strength located at the origin. 

The Laplacian growth equation Coming back to the derivation, we apply the conformal 
map from the unit disc in the complex plane w = exp(— p + i(f)), where the (stream) 
function (f){x,y) is harmonically conjugate to p{x,y), into the domain D(t) in the 
"physical" complex plane z = x + iy, and zero maps to zero with a positive coefficient. 
Denoting the moving boundary as z{t,l), where / is the arclength along the interface, 
one obtains 

K = Mztzi) = -dnP = di<P, (102) 



It is trivial to see that the chain of three equalities in (102) represent respectively the 



definition of Vn in terms of a moving complex boundary, z{t,l) (the first equality). 



the kinematic identity expressed by the last equation in the system (101) (the second 
one), and the Cauchy-Riemann relation (the last one) between p and (p. After 
reparamerization, / ^ we arrive to the equation 

lm{ztz^) = 1. (103) 



which possesses many remarkable properties, as will be seen below. The equation (103) 



is usually referred as the Laplacian growth equation (LGE) or the Polubarinova-Galin 



equation. In [Ml EH] it was noticed a fully unexpected feature of the equation (103): 
the boundary, z{t,(f)), taken initially as a polynomial of w = exp(i0), will remain a 
polynomial of the same degree with time-dependent coefficients during the course of 
evolution, so new degrees of freedom, describing the moving boundary, will not appear. 
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An even more remarkable observation concerning the equation (103), belongs to 



Kufarev 

to w = 



who found that a boundary taken as a rational function with respect 
exp(i0) will stay as such during the evolution. Moreover, he managed to 
integrate this dynamical system explicitly, and found first integrals of motion associated 
with moving poles and residues of the conformal map, z{t,exp[i(j))), describing the 
boundary. The authors [Ml EH |66] have however noticed that all the solutions obtained 
are short-lived, both because of instability and due to the finite volume of D(t), which is 
destined to shrink, because of a sink(s) located inside. We will address these interesting 
observations in detail in the second section of this chapter. 



LGE in the evolutionary form It is of help to present (103) in the evolutionary form 



defined as the dynamical system, where the time derivative constitutes the LHS and 



does enter the RHS. For this purpose we rewrite (103) as 



ztz^ = i + to, 

where to is real. Dividing both sides by l-z^P, we will obtain 

Zt i + to 



Taking conjugate form both sides and multiplying by i, we will have 



.Zt_ 



l + itp 



The LHS is the analytic function outside the unit disk in the w-plane. In accordance 
with the last equation, the real part of this analytic finction along the unit circumference 
equals |^,/,|~^. To recover the analytic function from the boundary value of its real part 
at the unit circle is a well-known procedure involving either the Hilbert transform or 
the Schwarz integral. The result is 

P 2^' 



IZt 



(104) 



,0 e*^ - (1 + e)e*<^ 

where an infinitesimally small positive e indicates correct limiting value of the integral 
while approaching the unit circumference. This useful formula was obtained by 
Shraiman and Bensimon in 1984 [67J. This expression for (103) in the evolutionary 



form reveals the nonlocal nature of Laplacian growth due to the integral in the RHS. 



The equation (104) helps to prove a beautiful statement, that every singularity of 



the function z{t, w) moves toward the unit circle from inside, or in other words the radial 
component of the 2D velocity of any singularity of the conformal map is positive. To 
prove the claim, we replace in (104) by W, defined earlier as Vl^ = —p + i(f). Then 
after we notice that 



zt{t,e 
'zwit.e^) 



'dW 



z=const 
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Figure 6. A DLA cluster, n = 100000. 



and that near a singular point w = a we can replace W = \og{w) = log(a), we can 
rewrite the real part of (104) as 
d log \a\ 



dt 



> 0. 



(105) 



Thus, we proved that each singular point of the conformal map moves toward the unit 
circle from inside, so the origin is a repellor for this dynamical system, and the unit 
circumference is an attractor. 



Diffusion limited aggregation The physics section of the survey cannot be completed 
without mentioning a fascinating discovery by T.A. Witten and L.M. Sander, who 
observed [58] in 1981 that a cluster on a 2D square lattice, grown by subsequent 
attaching to it a Brownian diffusive particles, eventually becomes a self-similar fractal 
(see Figure |6]) with a robust universal fractal (Haussdorff ) dimension given by 

Do = hm 1°^^ = 1.71 ± .01, (106) 
e->o log(l/e) 

where 1/e is a linear size of a cluster measured by a small "yard stick", e, and A^(e) is 
a minimal number of (small) boxes with a side e, which covers the cluster. 

Remarkably, this fractal appeared to be self-similar after appropriate statistical 
averaging. This means that its higher multi- fractal dimensions, Dg, defined as 

^-0 log(l/e) 

where Pi stands for a portion of a tiny box of a size of e, covered by the cluster under 
consideration, appear to be equal to each other, and to Dq, which is 1.71, as indicated 
above. Later, these findings were significantly clarified and refined in many respects, 
but the major challenge: how to calculate the universal dimension defined above still is 
an open question (see a relatively recent review [68j and references therein). 
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This problem is tightly connected to the Laplacian growth. Until very recently there 
were numerous claims that the DLA process is drastically different from the Laplacian 
growth, and even statements appeared that the DLA and fractals grown in Laplacian 
growth belong to different universality classes ^\ . However, the recent experiments by 
Fraud and Swinney [70] made crystal clear that the multi-fractal spectrum of a cluster 
grown in a viscous fingering process in a Hele-Shaw cell (that is a Laplacian growth) 
coincides with the DLA spectrum up to the margin accuracy of 1%, which is the maximal 
accuracy in these measurements. Thus, despite of its discrete and a stochastic nature, 



DLA can be understood by a continuum and deterministic Laplacian growth (101). 



Related problems Below is a list of physical problems connected with Laplacian growth. 

First, there is the so-called "singular" Laplacian growth, where a growing domain 
consists of needles with zero areas and divergent curvature at the tip. The mathematical 
description for this dynamics should be reformulated, since the gradient of pressure p 
diverges near moving needle tips, so the boundary velocity should be replaced by an 
appropriately regularized law. Interesting works by Derrida and Hakim [71], and by 
Feterson [72] in this direction deserve special attention. 

There is also a considerable amount of works in so-called nonlinear mean-field 
dynamics, where a phase field is involved, which gradually changes from unity in one 
of moving phases toward zero inside the second one [731 El] • Many of these processes, 
including dynamics of miscrostructure [75] in materials, growth of bacterial colonies 
in nutritional environment [76], and spinodal decomposition [77], governed by the 
time-dependent Ginzburg-Landau and the Cahn-Hilliard equations, are reduced to the 
Laplacian growth interface dynamics in a special singular limit, when the phase field 
degenerates to a step-function, thus becoming a characteristic function of a moving 
domain with a well defined boundary [751 179] . This is certainly worth to mention, both 
because it significantly enriches a physical process by introducing an additional field 
(the phase field) and since this is conceptually related to a random matrix approach 
to Laplacian growth, addressed in the survey, and where a distribution of eigenvalue 
support will play a role of a mean-field phase, introduced in this paragraph. 

Let us also mention several more "selection puzzles" , which belong to the Laplacian 
growth in various settings: selection of a shape of a separated inviscid bubble, observed 
by Taylor and Saffman in a viscous flow in a rectangular Hele-Shaw cell [80] from a 
continuous family of possible solutions (not to be confused with the Saffman- Taylor 
flngers family described in [51]); selection of a so-called "skinny" finger in a Hele-Shaw 
cell accelerated by a tiny inviscid bubble near the nose of a finger [81] ; and prediction of 
the periodicity for the so-called side-branching structure in dendritic growth [82j. These 
phenomena have the same (or almost the same) mathematical description. 

Another important comment about physics of Laplacian growth is that Darcy's 



law (92) is invalid near walls of a Hele-Shaw cell, including proximity to both parallel 
plates. This is because averaging of the Stokes fiow, yuV^v = Vp, given by (99) will no 



longer bring us to (100), due to boundary layer effects. This apparent difficulty gives 
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rise to tlie study of an interface dynamics with a Stokes flow, which is an extension of 
the Hele-Shaw (Darcy's) flow. The Stokes flow also contains remarkable physics and 
beautiful mathematics [831 EH ES] , which is still yet to be fully understood. 



4-3. Exact solutions 



(108) 

is conformal inside the 



Cardioid Consider the equation of motion for the droplet boundary under Laplacian 
growth 

lm{ztz^) = Q, 

where 27rQ stands for a rate of a source (sink). Here, z{t,e'^'^ 
unit circle, \w\ < 1, 
the solution, 

z = r{t)e''^, (109) 

comes to mind first, as the simplest one. It describes initially circular droplet centered 
at the origin, which uniformly grows (shrinks) while continuing to be a circle. Indeed, 



0^0, and w = e**^ in the equation. When one tries to solve (108) 



substituting (109) into (108) one obtains 



r(t) = y/2i\Q\T + Qt), (110) 

where a constant of integration, T stands for an initial time. When Q < (suction), 
the circle shrinks to a point at t = T, and the solution (109) ceases to exist after T. 
Could one find any other exact solutions, less trivial than given by (109)? 



Remarkably, the answer is yes, despite of nonlinearity of the Laplacian growth 



equation (108). Let's add to (109) an initially small quadratic correction. 



z = r{t)e'^ + a{t)e 



2i(j> 



111) 



The domain bounded by the curve described by (111), named a cardiod, is connected if 
\a\ < r/2. Substituting (111) into (108) one obtains two coupled nonlinear first order 
ODEs w.r.t. r and a: 

rr + 2aa = Q 
ar + 2 ar = 0, 



(112) 



with an easily found solution 

ar^ = ao 
r2 + 2a2 



2{\Q\to + Qt) ^^^^^ 

with ao and to as constants of integration. If Q > (injection), the cardiod will grow 
becoming m ore a nd more like a circle during the evolution. If instead Q < (suction) 
the cardioid ( 111 ) shrinks, deforms and ceases to exist after t* = to+3aJ^ /{^Q). This 



happens when the critical point of the conformal map given by (111) reaches the unit 



circle from outside. Then the cardioid ceases to be analytic and earns a needle-like cusp 
(a point of return with infinite curvature). This cusp is called type 3/2 (alternatively 
(2,3)-cusp) because in local Cartesian coordinates it is described by the equation y"^ ~ x^. 
We will see later that this kind of cusps is typical for those solutions of Laplacian growth 
which cease to exist in finite time. 
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Polynomials As a generalization, we are going to prove now that all polynomials of 
w, which describe boundaries of analytic domains when \w\ = 1, are solutions of (108). 
Assume a droplet is initially described by a trigonometric polynomial (with all critical 
points lying outside the unit disk, because its interior conformally maps onto a droplet): 



N 

E 

k=l 



ik<l> 



(114) 



Substituting (|TT4|) into (|T08|), one obtains coupled ODE's for time-dependent 
coefficients a^, and remarkably there are no other degrees of freedom which appear 
during the evolution. In other words, the evolving droplet will continue to be described 
by the polynomial (114), with coefficients, a^, changing in time in accordance with these 
ODE's: 

N-n 

[k ak hk+n + (fc + n)dfc ak+n] = QSn,o n = 0,1, . . . , N - 1. (115) 

k=l 



Moreover, ( |115[ ) can be integrated explicitly. Indeed, we notice first that the equation 
for k 



N — 1, namely 

ai (In + A^di = 0, 
is trivially solved with the answer 

Cn, 



N 

(l]\f 0,1 



;ii6) 



;ii7) 



where Cn is the constant of integration. Substituting (117) into the (A^ — 2)"^^ equation, 
which has a form 



ai a^r_i + (A^ — l)di ajv„i + a2 Aat + Na2 = 0, 



118) 



we notice that the LHS of (118) is proportional to a full derivative from the expression 



N- 



and is zero in accordance with the RHS of (118). Thus we obtain 



cln-i -\- N Cm a2 — Cn-i a\. 



(119) 



where C^-i is a constant of integration. Knowing a^^i and a^v in terms of ai and a2 
we can easily integrate the third equation from the end of the system (115), namely the 
(A^ — "iy^ equation. The result is 



aN-2 + {N-l)CN-ia2ai+NCNa^ai-N{N+l) 



Cn(^2 



C^_2at(120) 



Continuing in this way, we obtain an explicit dependence of linear 
combination of constants of motion, Ck, with coefficients which are polynomial forms 
w.r.t. ai, 02, . . .. The equation for n = from ( 115 ) already constitutes the full derivative 
and, as such, is trivially integrated: 

N 

Y,kW? = '^{Co + Qt), (121) 
fc=i 
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where Cq is a constant of integration. Here the LHS is a (scaled) area of the droplet, and 
the equation states that the area changes linearly in time. In other words, we integrated 
the system (115), and the solutions are polynomial forms with respect to a^, linear w.r.t. 
integrals of motion, C^, explicitly obtained. 

As in the case of cardioid, in case Q > the dynamics is stable and the droplet 
becomes eventually more and more round since all decay in time, but Oi in contrary, 
grows, as one can easily verify by looking to the system ( 115 ). If Q < 0, then the droplet 
shrinks and the solution ceases to exist in finite time. This happens because a critical 
point (s) hits a unit circle from outside manifesting a break of analyticity by making a 
cusp (of a 3/2 kind in general case). Except such circle centered at the 

location of sink, the solution stops to exist prior to the formation of a cusp, because of 
the droplet being completely sucked by the sink. 

The fact that a finite time singularity (a cusp) is generic follows directly from (117): 
since the conformal radius, Oi, should decrease as the area shrinks, then the coefficient. 
Oat, grows in time by virtue of (117), eventually bringing the system to a cusp. 

Now consider the external Laplacian growth, where an inviscid bubble, surrounded 
by a viscous fluid grows (shrinks) because of a source (sink) at infinity. Then we map 
conformally the exterior of the unit disk in the w-plane to the exterior of a bubble 
(viscous region) in the physical 2;-plane with a simple pole and positive residue (which 
is a conformal radius) at infinity. 



Here an analogy of the polynomial ansatz (115) will be the formula 



N 



-ik(f> 



(122) 



where a_i = r is the conformal radius, that is the radius of a circle perturbed by the 
rest of Ofc's. This case is also integrable in a way, very similar to the interior case shown 
above [HS]. One can also see that for an unstable LG, that is a growing bubble in the 
exterior problem, a finite time cusp is unavoidable. Indeed, the system of ODEs for the 
ansatz (122) will look the same as (115), but with values of n extended from — 1 to A^ + 1. 
Thus one can easily see that = Cn^^'^- This means that the highest harmonic will 
grow faster than a conformal radius, which should eventually break domain's analyticity 
through a cusp [87] . 

The area, to, of the growing bubble in this case equals 

t^ = \r\''-Y,k\ak\\ (123) 

fc>0 



Consider the simplest non-trivial example for (122), which describes a shape with 
three-fold symmetry. 



z{w) = rw + 



(124) 



we have 



3Cr2 



(125) 
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where C is a constant of integration, and the scaled area of the droplet identified with 
time t in this case, is: 

= t = \r\^ - 18\C\W\ (126) 

Clearly, this polynomial in |rp has a global maximum at solving 36|rcp|Cp = 1. 
We call the corresponding value of the area tc = j1 critical^ and conclude that the 
dynamics will lead to finite-time singularities for any initial condition ts 7^ Figure [7} 




Figure 7. A Hele-Shaw droplet approaching the critical area. 



In summary, we have shown that the Laplacian growth is integrable for polynomial 
(time-dependent) conformal mappings both in interior and exterior problem, and in 
both cases (growing of a bubble in the exterior and suction of a droplet for the interior 
problem) finite time singularities in the form of the 3/2-cusps are unavoidable. This is 
caused by an ill-posedness of the Laplacian growth without regularized factors, such as 
surface tension. 



Rational functions Kufarev [66] found a class of rational solutions of the equation ( 108 ) 
with simple poles. Let us show that all rational conformal maps from exterior (interior) 
of the unit disk to the exterior (interior) of a domain D are solutions of the LGE (108). 
We will include in our proof multiple poles for the sake of generality. Specifically, we 
claim that the expression 



Pk 

fc=l 1=1 



A 



kl 



ak 



\Pki 



(127) 



where pki G N, solves (108). Indeed, after substitution of (127) into (108), putting 
w = e^'^, we obtain a double sum, which we can decompose to elementary fractions with 



respect to (e**^ — a^)^*' by using repeatedly the identity 



(e^'^-afc)(e-^'^-a/) l-akOi 



1 + 



Di4> 



ai 



ak 



-i4> 



ai 



Equating coefficients prior to all independent modes to zero and sum of all constants 
(the zeroth mode) to Q in accordance with (108), we see, after some algebra, that all the 
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expressions are full derivatives, and after integration we obtain the following equations: 

- Eti E£i A,i[-zil/-a,)](^'^'~'y{p,i - 1)\ = C + Qt 



z{l/ak) = (3k 



k = 1,2, 



.N 



(129) 



A 



kl 



1,2, 



Oikl[ZwO-lw)Y'"\w=ak ' 

where C, an, and Pk are constants of integration. It is possible to show that in unstable 
LG, that is a an exterior problem with growth or an interior problem with shrinking, all 



the solutions (127) blow up in finite time by forming cusps, generally of the 3/2 kind. 



Another interesting class of rational solutions was also found by Kufarev [66j in case 
there are several sources instead of one, located at with rates Qk, and A; = 1, 2, . . . , A^. 
In this case the velocity potential (scaled pressure) diverges near Zk logarithmically with 
coefficients Qk- 

— p = Qk ^og{z — Zk) + regular terms (when z — > Zk)- (130) 
In this case the Laplacian growth equation has a form 

Qk 



lm{ztz^) = Re ^ — 



k=l 



bk{t)e 



(131) 



where bk{t) are time dependent inverse conformal pre-images of sources locations, Zk, so 
that 



M. 



Zk = z{h]^ 

In this case, the most general rational solution has a form 

A' Pk ^ M 

z = rw + -. r — h 



Bk 



k=l 



= 1 



k=l 



w - bk 



(132) 



(133) 



The result of integration is then given by (129), where summations incorporate the last 



sum in (133), the equation (132), and 

Bk = Ckt[z^a,il/w)]\^=i^ 



k = l,2. 



M, 



(134) 



where Ck are additional constants of motion. It is worth to mention that even if the 
initial configuration ;z(0,e*'^) does not include poles at bk (knows nothing about sources 
Qk at Zk), the solution earns terms with simple poles at bk immediately from the start, 
as one can see from the last equation. 

Thus, the singularities of any solution can be split into those imposed by the source 
location {bk in our case) and those determined by initial configuration, that are a^. 

One should also beware that the interface can reach sources during evolution, thus 
breaking analyticity by forming a cusp, and after this moment a solution ceases to exist. 



Logarithms In the paper [88] Kufarev and Vinogradov have found a logarithmic class 



of solutions of (103), which was later rediscovered and studied in detail by several 



authors [89l [90], IHH [921 193]. This class appeared to be particularly fruitful from both 
mathematical and physical points of view: besides providing a significant extension from 
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rational solutions, the logarithmic ones are often free of finite time singularities for an 
unstable exterior problem, which is the most important for physics. Existence of these 
solutions for all times allows to study the long time asymptotics, which is perhaps the 
major goal of this research. These so-called multi-logarithmic solutions have a form 

N 



rw 



( 

+ ^ afc log [ 

k=i ^ 



w 



- 1 



(135) 



where ra^, are parameters (some of them are time dependent), and \ak\ < 1 for a 
conformal mapping from the exterior of the unit disk. Using the method outlined above 
for rational solutions one could easily figure out that (135) satisfies the LGE (103) with 
all ttfc to be constants in time and the following time dependence of r and a^: 



'"^ + EfeLi Ez=i "fc log(l - ttkai) = C + Qt 



(136) 

r/ak + E£i log(;4 -l)=Pk; k = l,2,...,N, 

where C and f3k are constants of motion. It is less trivial to show that the solutions 
(135) may be free of finite time singularities, but the following example illustrates it well: 

a exp {27iik/N) and 



let's impose a Zat symmetry over the system (135) by setting ak 

= a exp {2nik/N), with positive a and a. Then (136) looks significantly simpler 



r2 + iV«2 ^^v^^ iQg^i _a^^,) = C + Qt 
r/a + a J2k=i 7fc log(^ - 1) = /5, 



(137) 



where = e^'^^^^^ . Equating the derivative of (135) to zero, we find critical points, hk- 
As expected, bk = b'jk and 



(138) 



Assuming the initial b to be positive, we see that b is always positive, if a is, which is 
always the case, since as follows from the second equation in (137) 

^2N-2 



- + 2aN- 
a 1 



.2N 



(139) 



and therefore a is positive, since r is. Let's also notice that a cannot reach 1 since this 
would make the RHS of the second equation in (137) infinite, which would contradict 
to the fact that it is a finite constant. Thus, from 

< 6 < a < 1, d > 0, and b>0 (140) 



it follows that critical points and singularities of the conformal map (135) will always 
stay inside the unit circle, which guarantees existence of the solution (135) for all times. 
This simple example illustrates the fact that many of these solutions are free of finite 
time-singularities (see details in [^21122]), but the interesting problem of comprehensive 



classification of initial data for the solutions (135) which do not blow up in finite time 
still is an open question. 
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4.4- Mathematical structure of Laplacian growth 

4.4-1- Conservation of harmonic moments The following remarkable property of the 
Laplacian growth was found by S. Richardson in 1972 [91] for a point-like source Q at 
the origin, for which 

He showed that all positive harmonic moments of the viscous domain, -D(t), 



/ dxdy, k = 1,2, 

JD(t) 



I ,kVn-= i ipd^z' - z^ dnP) -, (142) 

JdD(t) ^ JdD(t) ^ 



Ck = I z^dxdy, A; = 1,2,... (141) 

m 

do not change in time, while the zeroth moment, which is the area of the growing bubble, 
changes linearly in time: 

dCk f i.^. dl f , ^ u I. ^ . dl 

dt JdD(t) " JdD{t) 

[dl is an element of arclength) because Vn = —dnP and p = along the boundary dD(t). 
By virtue of Gauss' theorem, it equals 

/ WipVz" - z'Wp) ^ = Q 4 0. (143) 

This property may be used as the definition of the idealized Laplacian growth problem, 
namely to find an evolution of the domain whose area increases in time, while all positive 
harmonic moments do not change. 

4.4- 2- LG and the Inverse Potential Problem One can easily notice that the harmonic 
moments are the coefficients of the (negative) power expansion of the so called Cauchy 
transform Cr){z) of the domain D, namely 

Co(z) = i / = f ("^' 

TT In z — z ^-^ z*^^^ 

Since the Cauchy transform, Cd{z), is the derivative of the Newtonian potential $(z) 
created by matter occupied the domain D with a unit density, 

$(z) = I log|z-z'|^^, (145) 
Jd 

we see a deep connection between the Laplacian growth with the so-called inverse 
potential problem, asking to find a domain D occupied uniformly by matter which 
produces a given far field Newtonian potential. The harmonic moments in this context 
are multipole moments of this potential. If the domain D{t) grows in accordance with 
the idealized Laplacian growth, then the potential ^{z) changes linearly in time, so (up 
to a constant): 

$(t,^) = |Jlog|z|, (146) 
which is a potential of a point-like (increasing in time) mass at the origin. 
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4.4-3. Laplacian growth in terms of the Schwarz function Let F{x,y) = define an 
analytic closed Jordanian contour F on the plane. Replacing the cartesian coordinates, 
X and y, by complex ones, z = x + iy and z = x — iy, one obtains a description of F as 

z + z z — z\ ^, ^ , ^ 

Solving the last equation with respect to z one obtains: 

z = S{z), (148) 

when 2 G F. The function S{z) is called the Schwarz function of the curve F [05]. It is 
the same mathematical object we encountered in the previous chapter. This function 
plays an outstanding role in the theory of quadrature domains (see next chapter). It 
has the following Laurent expansion, valid at least in a strip around the curve F: 

^(-) = E ^ + E (149) 

k=0 k=0 

where tk are the external harmonic moments defined as 

U-\f k = l,2,..., (150) 

where D_ is the domain complimentary to the domain D. 



From (144) and (149) we obtain the connection between the Cauchy transform of 



a domain with the Schwarz function of its boundary: 

/ S(z') dz' 

Rewriting the Laplacian growth dynamics in terms of the Schwarz function, S{z) 
[87] one obtains 

dtS{z) = 2dM. (152) 

where W = —p+i(j) = log w is the complex potential defined earlier. This last form of the 
Laplacian growth is very instructive. In particular, it helps to understand the origin of 
constants of integration in all exact solutions of the Laplacian growth equation presented 
above as a result of direct integrating efforts. Indeed, the RHS in the last equation is 
analytic in the viscous domain D(t) except a simple pole at the origin (we consider 
an internal LG problem with a source at the origin). In order for the LHS to satisfy 
this condition, all the singularities of S{t, z) outside the interface should be constants of 
motion. At zero the Schwarz function should have a simple pole with a residue (which is 
the area of the domain D{t)) linearly changing in time. This observation can be easily 
seen as an alternative proof of the Richardson theorem, stated above. 
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4-4-4- The correspondence of singularities The Schwarz function is connected to a 
conformal map z = /{w) from the unit circle to the domain D through the following 
formula [95] 



S{z) = f 



(153) 



where f~^{z) is the inverse of the conformal map w = f~^{z). This formula helps to 
derive a one-to-one correspondence between singularities of S{z) inside D and f{w) 
inside the unit circle: if near a singular point a the conformal map f{w) diverges as 

A 



/H 



[w 



(154) 



(here by convention p = stands for a logarithmic divergence), then the Schwarz 
function S{z) diverges near a point b = /(1/a) with the same power, p, as 

B 



where 



Siz) 



A 



by 



B 



(155) 



(156) 



w=l/a 



B and b are constants of motions as showed above, thus the last formula together 
with the relation b = /(1/a) and the area linearly changing in time and expressed in 
terms of the parameters of fiw) constitute the whole time dynamics of singularities of 
f(t,w) [SI [96|. The reader can see the equivalence of these formulae with constants 
of integration obtained earlier when various classes of exact solutions were derived by 
direct integration. 



4.4-5. A first classification of singularities As mentioned in the previous sections, 
existence of the singular limit was established at the same time with the model [Ml ES] . 
It became a fertile field of study in itself, and led to further developments of the problem 
[SH 197]. In a series of papers [HHl [991 HOH HDSl IDS], the possible boundary 
singularities were studied, as well as the problem of continuing the solutions for certain 
classes. It was found that, in the free-space set-up, the generic critical boundary features 
a cusp at {xQ,yQ), with local geometry of the type 

(x — XqY ~ (y — VoY, {p, q) mutual primes. (157) 

The most common cusp is characterized by g = 2,p = 3, but g = 2,p = 5 can also 
be obtained fairly easy by choosing proper initial conditions. Very special situations, 
where a finite-angle geometry is assumed as initial condition were also considered |1U4] . 

It was shown be several methods that dynamics can be continued through a cusp 
of type (2, 4A; + 1), A; > [Ml HUS]. 
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4.4-6. Hydrodynamics of LG and the singularities of Schwarz function As indicated 
above, the Schwarz function encodes information about the conserved moments {tk}, 
through its expansion at infinity [50] : 

to 

z 



S{z) 



J2tkz'-' + 0{z-'). 



(158) 



k>0 



This function is useful when computing averages of integrable analytic functions f{z) 
over the domain (an interior domain): 

M „ 

f{z)dxdy = > ] y]c,kf^'^{zk) + V / hm{z)f{z)dz, (159) 



k=l i=l 



+ ^ / hm{z)f{z)dz, 
Zk and branch cuts 7„ 



if the function 5* has poles of order n/^ at z 
functions hm{- 

f{z) = Xd+{z) and taking a derivative with respect to t^, we obtain the relation 



with jump 



Applying formula (159) for the characteristic function of the domain 



d 

dto 



y^^ResSjzk) + / hjn{z)dz 



1, 



(160) 



which shows that the singularity data of the Schwarz function in can be interpreted 
as giving the location and strength of fluid sources (isolated or line-distributed) [M] . 
Identifying the 2D uniform measure with another, singular (point or line-distributed) 
distribution, is referred to as sweeping of a measure. We will repeatedly encounter this 
process in the next chapter. In the case when the Schwarz function is meromorphic in 



(it has only isolated poles as singular points), ( |159[ ) becomes 
f{z)dxdy 



N rife 



(161) 

k=i i=i 

and the domain is called a quadrature domain |106l I107|. 11081 1109i 1110] . Generically, the 
Schwarz function may have branch cuts in D^, in which case is called a generalized 
quadrature domain [lllj. This is the typical scenario for our problem. The rigorous 
theory of quadrature domains is outlined in the next chapter. 



The hydrodynamic interpretation of the Schwarz function arises from(152), which 
is worthwhile to rewrite here 



dtSiz)= dM 



(162) 



after rescaling by 2. Let C be some closed contour, boundary of a domain B, and 
integrate equation ( 1162 ) over it. We obtain 



dt (t S{z)dz = // uj dxdy — i jj Vvdxdy, (163) 
Jc JJ B JJb 

where u = dyV^ — d^Vy is the vorticity field, and Vv = d^v^ + dyVy is the divergence of 
velocity field. The real part of this identity shows if the fiow has zero vorticity, we have 



Re dt j S{z)dz = 0. 



(164) 
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The imaginary part of (163) illustrates again the interpretation of singularity set of S{z) 



as sources of water (which occupies Z)+ in a canonical Laplacian growth formulation, 
while the exterior domain, D_, is occupied by a viscous fluid, which we call oil [50]): 



assume that the contour C in (163) encircles the droplet without crossing any other 



branch cuts, then the contour integral may be performed using Cauchy's theorem, giving 
the total flux of water: 



V ■ V dxdy = Q = 1. 



We note here that equation ( 162 ) implies existence of a closed form 
dVt = Sdz + Wdt, 
whose primitive Q has for real part the Baiocchi transform of p: 



Re n 



p{z, T)dT. 



(165) 



(166) 



(167) 







One can see that Re Q coincides with the potential $ introduced earlier. From the 
continuity equation for water p + Vv = and the Darcy law for water (opposite to oil) 
V = Vp, we obtain for the time evolution of water density at a given point z, 

p = -Ap^ p{z,t) = p{z,0) - A [ p{z,T)dT. (168) 



Equation (168) may be immediately generalized in a weak sense, replacing the water 
density by the characteristic function of the domain -D+,p — > which shows that 
the Baiocchi transform Re Q may be interpreted as the electrostatic potential giving 
the growth of the water domain. 



Similarly, applying an antiholomorphic derivative to (162), we obtain 

-Ap + iA(j), 



V ■ V 



(169) 



so that the imaginary part of the form can be considered an electrostatic potential 
for the time integral of vorticity at a given point z: 

AImfi(2,t)= / Lj{z,T)dT. (170) 



4-4-7- Variational formulation of Hele-Shaw dynamics Formula (159) has another 



physical interpretation, which we explore in this section. Besides hydrodynamics, it 
also allows to describe the droplet through a variational (minimization) formulation, 
which will become very relevant when considering the singular limit. 

Consider the case when the Schwarz function has only simple poles {zk} and cuts 
at {'Jm}, with residues Res S{zk) and jump functions hm{z), inside the droplet. A 
simple calculation shows that these singular points constitute electrostatic sources for 
the potential Re Q: 



A Re n{z) = Res S{zk)5{z - Zk) + ^ 



,{C)5{z - QdC. (171) 
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If we apply (159) to all positive powers z^,k > 0, we conclude that the singular 
distribution {-Zfc}, {7m} and the uniform distribution p{z) = Xd+{z) have the same 
interior harmonic moments Vk = {z'^),k > 0. Thus, they create the same electrostatic 
potential outside the droplet. It is therefore possible to substitute the actual singular 
distribution {2fe},{7m} with the smooth, uniform distribution p{z) in calculations 
related to the exterior potential. Beyond the mathematical equivalence, however, this 
fact has an important physical interpretation, whose full meaning will become apparent 
in the critical limit: when one more quantum of water is pumped into the droplet, it first 
appears as a new singular point of the Schwarz function (a 5-function singularity). After 
a certain time, though, the droplet adjusts to the new area (subject to the constraints 
given by the fixed exterior harmonic moments), and reaches its new shape (with uniform 
density of water inside). Therefore, we can say that the singular distribution {zk}, {7m} 
represents the fast-time distribution of sources of water, while the uniform distribution 
p{z) is the long-time, equilibrium distribution of the same amount of water. When the 
dynamics becomes fully non-equilibrium (after the cusp formation), this equivalence 
breaks down, and the correct distribution to work with is the set of poles (cuts) of 
the Schwarz function. In that case, the issue becomes solving the Poisson problem 
ARei7 = ResS'(2;fc)(5(-2 — -Zfc), and finding the actual (time-dependent) location of the 
distribution of charges Zk(t), subject to usual conditions for the electrostatic potential 

In the equilibrium case, however, it is appropriate to work with the smooth 
distribution p{z). Since the actual electrostatic potential ARef2 contains the regular 
expansion V{z) = Ylik^kZ^^ ^^^^ it to the contribution due to the distribution 
p{z). We obtain for the total potential: 

Hz,z)= [ p{C)\og\z-C\'d'C + V{z) + Vi7). (172) 
Jd+ 

Inside the droplet, this potential solves the Poisson problem A ^{z, z) = p{z) = 1, and 
on the boundary it creates the electric field E{z) = = z. This means that inside the 
droplet, this potential is actually equal to Therefore, the problem of finding the 
actual shape of a droplet of area to and harmonic moments {t^} can be stated as: 

Find the domain of area to such that ^{z, z) = l^p on 

Since p{z) is the characteristic function of D^, we may also write this problem in the 
variational form: 



5p{z) 



[ p(z)\\z\'-Viz)-Viz)- [ p(C)log|^-C|VC 
Jd+ L Jd+ 



d^z = 0. 



This equation is simply the minimization condition for the total energy of a distribution 
of charges p{z), in the external potential W{z,z) = — + V{z) + V{z). Therefore, 
the equilibrium (long-time limit) distribution of water has the usual interpretation of 
minimizing the total electrostatic energy of the system. However, when the system is 
not in equilibrium, this criterion cannot be used to select the solution. 
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5. Quadrature Domains 

We have seen in the previous sections that polynomial or rational conformal mappings 
from the disk have as images planar domains which are relevant for the Laplacian 
growth (with finitely many sources). The domains in question were previously and 
independently studied by mathematicians, for at least two separate motivations. First 
they have appeared in the work of Aharonov and Shapiro, on extremal problems of 
univalent function theory |112] . About the same time, these domains have been isolated 
by Makoto Sakai in his potential theoretic work [113J. These domains, known today as 
quadrature domains, carry Gaussian type quadrature formulas which are valid for several 
classes of functions, like integrable analytic, harmonic, and sub-harmonic functions. 
The geometric structure of their boundary, qualitative properties of their boundary 
defining function, and dynamics under the Laplacian growth law are well understood. 
The reader can consult the recent collection of articles |114] and the survey |115] . The 
present section contains a general view of the theory of quadrature domains, with special 
emphasis of a matrix model realization of their defining function. 

This chapter is organized in the following way: after presenting the theory of 
quadrature domains for subharmonic and analytic functions, we give an overview of 
the (inverse) Markov problem of moments, followed by its analogue in two dimensions, 
which is based on the notion of exponential transform in the complex plane. The 
following sections illustrate the reconstruction algorithm for the shape of a droplet, and 
point to a few essential properties of the problem for signed measures. 

5.1. Quadrature domains for subharmonic functions 

Let be a subharmonic function defined on an open subset of the complex plane, that 
is Aip > 0, in the sense of distributions, or the submeanvalue property 



holds for any disc centered at a, of radius r, B{a, r) contained in the domain of definition 
of (p. Henceforth dA denotes Lebesgue planar measure. Thus, with Q = B{a,r), 
c = \B{a,r)\ = 7rr^ and fi = c6a there holds 



for all subharmonic functions (p in fi. This set of inequalities is encoded in the definition 
that f2 is a quadrature domain for subharmonic functions with respect to n |113j , and 
it expresses that fl = B{a,r) is a swept out version of the measure /i = cSa- If c 
increases the corresponding expansion of f2 is a simple example of Hele-Shaw evolution, 
or Laplacian growth, as we have seen in the previous section. 

The above can be repeated with finitely many points, i.e., with /x of the form 





(173) 



H = Ci6ai + . . . + CnS, 



(174) 



5.1 Quadrature domains for subharmonic functions 



47 



ttj G C, Cj > 0: there always exists a unique (up to nullsets) open set 1] C C such that 



(173) holds for all ip subharmonic and integrable in Q. One can think of it as the union 
Uj=i B{aj, rj), Tj = ^Jcj/Ti, with all multiple coverings smashed out to a singly covered 
set, VL. In particular, (J"^^ _B(aj, rj) C VL. 

The above sweeping process, i— >• f2, or better ^ ^ xn- (dA), called partial balayage 
[113] . |116j . |117j . applies to quite general measures /i > and can be defined in terms 
of a natural energy minimization: given fi, u = xn ■ (dA) will be the unique solution of 



Minimize!, \ \^ — v\\e s.t. v < dA, J ~ J ^f^- 
Here || ■ ||e is the energy norm: 

1 f 1 

= (/W,At)e, with (/i, i/)e = — / log^ -dfi{z)du{C). 



If fi has infinite energy, like in (174), one minimizes — 2(/x, u) + ||z/||g instead of — 
which can always be given a meaning [26] . 
By choosing 

^{0 = ±log\z - C\ 



2 



(175) 



in equation (173), the plus sign allowed for all z G C, the minus sign allowed only for 
z ^ fl, one gets the following statements for potentials: 

> in all C, 
Uf" = outside n. 

Here 

denotes the logarithmic potential of the measure /i, and = JJ^"''^^. In particular, 
the measures /i and xn ■ (dA) are gravi- equivalent outside Q. By an approximation 



argument, (175) is actually equivalent to (173). 



Let us consider now an integrable harmonic function h, defined in the domain 
Since both ip = ±h are subharmonic functions, we find 

J hdA = J hdfi = '^^Cjh{aj). (176) 

That is, a Gaussian type quadrature formula, with nodes {aj} and weights {cj} holds. 
We say in this case that is a quadrature domain for harmonic functions. Similarly, one 
defines a quadrature domain for complex analytic functions, and it is worth mentioning 
that the inclusions {QD for subharmonic functions} C {QD for harmonic functions} C 
{QD for analytic functions} are strict, see for details |113] . 

Recall that for a given positive measure a on the line, rapidly decreasing at infinity, 
the zeros of the A^-th orthogonal polynomial are the nodes of a Gauss quadrature 
formula, valid only for polynomials of degree 2A^ — 1. The difference above is that 
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the same finite quadrature formula is valid, in the plane, for an infinite dimensional 
space of functions. A common feature of the two scenarios, which will be clarified in 
the sequel, is the link between quadrature formulas (on the line or in the plane) and 
spectral decompositions (of Jacobi matrices, respectively hyponormal operators). 

Let K = convsupp/i be the convex hull of the support of fj,, i.e., the convex hull 
of the points ai, . . . , a„. As mentioned, fl can be thought of as smashed out version 
of IJj=i ^{^ji^j)- The geometry of Vt which this enforces is expressed in the following 



sharp result ( |116j . |118] . |119] ): assume that Vt satisfies (173) for a measure > of 



the form (174). Then: 



(i) dVL may have singular points (cusps, double points, isolated points), but they are 
all located inside K. Outside K, dfl is smooth algebraic. 

For z G dQ \ K, let denote the inward normal of dQ at z (well defined by (i)). 

(ii) For each z G dQ \ K, intersects K. 

(iii) For z,w E dQ \K, z ^ w, and A'^ do not intersect each other before they reach 
K. Thus Q\K is the disjoint union of the inward normals from dQ \ K. 

(iv) There exist r(z)>OfoTzEKr\Q such that 



n= \J B{z,r{z)). 



zeKnn 

(Statement (iv) is actually a consequence of (iii).) 

To better connect our discussion with the moving boundaries encountered in the 
first part of this survey, we add the following remarks. Since Q is uniquely determined 
by {aj,Cj) one can steer Q by changing the Cj (or aj). Such deformations are of 
Hele-Shaw type, as can be seen by the following computation, which applies in more 
general situations: Hele-Shaw evolution Q{t) corresponding to a point source at a G C 
("injection of fluid" at a) means that Q{t) changes by dfl{t) moving in the outward 
normal direction with speed 

<9Gn(t)(-,a) 



dn 

Here Gn{z,a) denotes the Green function of the domain Q. If (f is subharmonic in a 
neighborhood of Q(t) then, as a consequence of G^{-,a) > 0, Gfi{-,a) = on dQ and 
-AGn{-,a) = 6a, 

— f (pdA = [ (speed of dVL{t) in normal direction) (^ds 
dt Jn{t) Jdn{t) 

dGnit){-,a) f dip 



V9AGn(t)(-,a)dA+ / Gn(i)(-,a) Ay^dA > </?(a). 

Q.{t) Jn{t) 
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Hence, integrating from t = to an arbitrary t > 




> 



j 

Jn 



n{o) 



Lp dA + t{p{a), 



telling that if ^2(0) is a quadrature domain for fi then Q{t) is a quadrature domain for 
fx + Ua- 

We remark that quadrature domains for subharmonic functions can be defined in 
any number of variables, but then much less of their qualitative properties are known, 
see for instance jl 14] . 

5.2. Quadrature domains for analytic functions 

Critical for our study is the regularity and algebraicity of the boundary of quadrature 
domains for analytic functions. This was conjectured in the early works of Aharonov 
and Shapiro, and proved in full generality by Gustafsson |108j . A description of the 
possible singular points in the boundary of a quadrature domain was completed by 



Assume that the quadrature domain for analytic functions Vt has a sufficiently 
smooth boundary F. Let us consider the Cauchy transform of the area mass, uniformly 
distributed on f2: 



This is an analytic function on the complement of f2, which is continuous (due to the 
Lebesgue integrability of the kernel) on the whole complex plane. In addition, the 
quadrature identity implies 



Therefore, by standard arguments in function theory one proves that the continuous 
function w extends meromorphically from F to Vt. The poles of this meromorphic 
extension coincide with the quadrature nodes. 

The converse also holds, in virtue of Cauchy's formula: if / is an integrable analytic 
function in fi, then 



Thus, we recover the following fundamental observation: if Vt is a hounded planar 
domain with sufficiently smooth boundary T , then Q is a quadrature domain for analytic 
functions if and only if the function w ^ w extends meromorphically from T to Q. 



Sakai [98ll99l [UK)] . 




n 




From the Stokes formula. 
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Note that above, and elsewhere henceforth, we do not assume that the weights 
in the quadrature formula for analytic functions are positive. In this way we recover 
the fact (already noted in the previous chapters) that quadrature domains for analytic 
functions are characterized by a meromorphic Schwarz function, usually denoted S{z). 
A second departure from the quadrature domains for subharmonic functions is that the 
quadrature data {aj,Cj) do not determine the quadrature domain for analytic functions. 
Indeed, consider the annulus Ar^R = {z, r < \z\ < R}. Then 

/ fdA = n{R'-r')f{0), 

for all analytic, integrable functions / in A^^r. 

The question how weak the smoothness assumption on the boundary F can be to 
insure the use of the above arguments has a long history by itself, and we do not enter 
into its details. Simply the existence of the quadrature formula and the fact that the 
boundary is a mere continuum implies, via quite sophisticated techniques, the regularity 
of r. See for instance pifM)]. 

The Schwarz function is a central character in our story. It can also be related to the 
logarithmic potentials introduced in the previous subsection. More specifically, given 



any measure as in (174) and any open set f2 containing supp fi, define (as distributions 
in all C) 

oz 

Then 

Au = xn- fJ', = 1 - Xn + fJ'- 



Note that with /i of the form (174) w is harmonic in f2 except for poles at the points Uj 



and that in particular, S{z) is meromorphic in f2. 



It is clear from (175) that is a subharmonic quadrature domain for if and only 
if M > everywhere and u = outside Q. Then also Vm = outside Q. Similarly, 
the criterion for Q being a quadrature domain for harmonic functions is that merely 
u = Vm = on C \ (The vanishing of the gradient is a consequence of the vanishing 
of u, except at certain singular points on the boundary.) To be a quadrature domain 
for analytic functions it is enough that just the gradient vanishes, or better in the 
complex- valued case, that || = on C \ (or just on dVl). 

Gustafsson's innovative idea, to use the Schottky double of the domain, can be 
summarized as follows. Let f2 be a bounded quadrature domain for analytic functions, 
with boundary F. We consider a second copy Q of Q, endowed with the anti-conformal 
structure, and "glue" them into a compact Riemann surface 

X = nuTun. 

This (connected) Riemann surface carries two meromorphic functions: 
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Any pair of meromorphic functions on X is algebraically dependent, that is, there exists 
a polynomial Q{z^w) with the property Q{g, /) = 0, and in particular 

Q{z,S{z)) = Q{z,z) = {], zeV. 

The involution (flip from one side to its mirror symmetric) on X yields the Hermitian 
structure of Q\ 



Q{z,w) = 



CLijZ W"^ ^ (Xij Ojji- 



One also proves by elementary means of Riemann surface theory that Q is irreducible, 
and moreover, its leading part is controlled by the quadrature identity data: 

Q{z,z)-\P{z)\^ = 0{z''-\z''-^), 

where 

P{z) = (z - ai){z - a2)...{z - an). 



This Riemann surface is the continuum limit of the spectral curve (67), for oo. 
Following Gustafsson ( |108] ). we note a surprising result: 

a) The boundary of a quadrature domain for analytic functions is a real algebraic, 
irreducible curve. 

b) In every conformal class of finitely connected planar domains there exists a 
quadrature domain. 

c) Every bounded planar domain can be approximated in the Haudorff distance by 
a sequence of quadrature domains. 

The last two assertions are proven in Gustafsson's influential thesis |108] . Recently, 
considerable progress was made in the construction of multiply connected quadrature 
domains, see [TMimifm] . 

It is important to point out that not every domain bounded by an algebraic curve is 
an algebraic domain in the above sense. In general, if a domain f2 C C is bounded by an 
algebraic curve Q{z,z) = (Q a polynomial with Hermitian symmetry), then one can 
associate two compact symmetric Riemann surfaces to it: one is the Schottky double 
of fl and the other is the Riemann surface classically associated to the complex curve 
Q{z,w) = 0. For the latter the involution is given by {z,w) i-^ (w,z). In the case of 
algebraic domains (this is another circulating name for quadrature domains for analytic 
functions), and only in that case, the two Riemann surfaces canonically coincide: the 
lifting 

z ^ {z,S{z)) 

from Q to the locus of Q{z, w) = extends to the Schottky double of Q and then gives 
an isomorphism, respecting the symmetries, between the two Riemann surfaces. 
As a simple example, the Schottky double of the simply connected domain 



n = {z = x + iy eC:x^ + y^ <l} 
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has genus zero, while the Riemann surface associated to the curve x"^ + = 1 has genus 
3. Hence they cannot be identified, and in fact Q is not an algebraic domain. 

Other ways of characterizing algebraic domains, by means of rational embeddings 
into n dimensional projective space, are discussed in [122j. 

5.3. Markov's moment problem 

We pause for a while the main line of our story, to connect the described phenomenology 
with a classical, beautiful mathematical construct due to A. A. Markov, all gravitating 
around moment problems for bounded functions. 

The classical L-problem of moments (also known as Markov's moment problem) 
offers a good theoretical framework for reconstructing extremal measures n from their 
moments, or equivalently, from the germ at infinity of some of their integral transforms. 
The material below is classical and can be found in the monographs |123l I124j . We 
present only a simplified version of the abstract L-problem, well adapted to the main 
themes of this survey. 

Let K he a compact subset of M" with interior points and let A C N" be a 
finite subset of multi-indices. We are interested in the set of moment sequences 

«(t(/) = / x'^ f{x)dx, a E A, 

of all measurable functions / : K — > [0,1]- Regarded as a subset of M'^I, is a 
compact convex set. An — L°° duality argument (known as the abstract L-problem 
of moments) shows that every extremal point of Tja is a characteristic function of the 
form X{p<7}) where we denote: 

{p < 7} = {x G K] p{x) < 7}. 

Above 7 is a real constant and p is an A-polynomial with real coefficients, that is 
p{x) = Xlo-eA Co-x*^. Indeed, to find the special form of the extremal functions /, one has 
to analyze when the inequality 

/ P{x)f{x)dx < \\p\\i\\f Woo = / \p{x)\dx 

is an equality. For a complete proof the reader can consult Krein and Nudelman's 
monograph |124j . 

As a consequence, the above description of the extremal points in the moment set 
Ha implies the following remarkable uniqueness theorem due to Akhiezer and Krein: 

For each characteristic function x of a level set in K of an A-polynomial there exists 
exactly one class of functions f in L'^{K) satisfying a{f) = a{x). For a non-extremal 
point a{f) e there are infinitely many non- equivalent classes in L°°{K) having the 
same A-moments. 
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Let us consider a simple example: 

K={{x,y)- x2 + y2<l}cM^ 

and 

fi+ = {(x, y)eK- X > 0, y> 0}, = {{x,y) e K- x < 0, y< 0}. 

The reader can prove by elementary means that the sets VL± cannot be defined in the 
unit ball K hy a. single polynomial inequality. On the other hand, the set 

Vt = VL+UVL_ = {{x,y)] xy > 0}, 

is defined by a single equation of degree two. 

Thus, no matter how the finite set of indices A C is chosen, there is a continuum 
fs, s G M, of essentially distinct measurable functions fs : K — > [0, 1] possessing the 
same A-moments: 

/ x''^y''^fs{x,y)dxdy= I x^^y^^dxdy, seR,aeA. 
J K Jn+ 

On the contrary, if the set of indices A contains (0,0) and (1,1), then for every 
measurable function / : K — ^ [0, 1] satisfying 

/ x'''y''^f{x,y)dxdy= / x^'y^^dxdy, a G A, 

J K Jq 

we infer by Akhiezer and Krein's Theorem that f = Xo,-, almost everywhere. 

On a more theoretical side, we can interpret Akhiezer and Krein's Theorem in 
terms of geometric tomography, see |125] . Fix a unit vector uj G M", ||ci;|| = 1, and 
let us consider the parallel Radon transform of a function / : K — > [0, 1], along the 
direction a;: 

{Rf){uj,s)= I f{x)dx. 

J {x,U!)=S 

Accordingly, the k-th moment in the variable s of the Radon transform is, for a 
sufficiently large constant M: 

{Rf)iuj,s)s''ds = / {x,iu)^f{x)dx = (177) 



-M JK 



^ M! /■ x'^u;^f{x)dx = J2 ^^'^^.l/)- (178) 

\a\=k ^' \a\=k ^' 

Since there are N{n, d) = C"^^ linearly independent polynomials in n variables of 
degree less than or equal to d, a Vandermonde determinant argument shows, via the 
above formula, that the same number of different parallel projections of the "shade" 
function / : K — > [0, 1], determine, via a matrix inversion, all moments: 

aaif), \cr\ < d. 
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The converse also holds, by formula (178). These transformations are known and 



currently used in image processing, see for instance |126j and the references cited there 
In 

process 



In conclusion, Akhiezer and Krein's Theorem asserts then that in the measurement 



only black and white pictures, delimited by a single algebraic equation of degree less 
than or equal to d, can be exactly reconstructed. Even when these uniqueness conditions 
are met, the details of the reconstruction from moments are delicate. We shall see some 
examples in the next sections. 

5.3.1. Markov's extremal problem and the phase shift By going back to the source 
and dropping a few levels of generality, we recall Markov's original moment problem 
and some of its modern interpretations. Highly relevant for our " quatization" approach 
to moving boundaries of planar domains is the matrix interpretation we will describe 
for Markov's moment. Again, this material is well exposed in the monograph by Krein 
and Nudelman jl24] . 

Let us consider, for a fixed positive integer n, the L-moment problem on the line: 



ak = akif) = f t''f{t)dt, 0<k<2n, 



where the unknown function / is measurable, admits all moments up to degree 2n and 
satisfies: 

0< f <L, a.e.. 

As noted by Markov, the next formal series transform is quite useful for solving 
this question: 

1 / oq ai a2n 



exp 



bo bi 

1 + - + 4 + • • ■ ■ 179 



Remark that, although the series under the exponential is finite, the resulting one might 
be infinite. 

The following result is classical, see for instance |123j pp. 77-82. Its present form 
was refined by Akhiezer and Krein; partial similar attempts are due, among others, to 
Boas, Ghizzetti, Hausdorff, Kantorovich, Verblunsky and Widder, see |123l 1124] . 

(Markov) Let ao, ai, . . . , a2n be a sequence of real numbers and let bo, bi, . . . be its 
exponential L-transform. Then there is an integrable function /, < / < L, possessing 
the moments ak{f) = ak, < A; < 2n, if and only if the Hankel matrix {bk+i)'^i^Q is non- 
negative definite. Moreover, the solution f is unique if and only if det{bk+i)k i=o = 0- -^''^ 
this case the function f /L is the characteristic function of a union of at most n bounded 
intervals. 

The reader will recognize above a concrete validation of the abstract moment 
problem discussed in the previous section. 
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In order to better understand the nature of the L-problem, we interpret below 
the exponential transform from two different and complementary points of view. For 
simplicity wc take the constant L to be equal to 1 and consider only compactly supported 
originals /, due to the fact that the extremal solutions have anyway compact support. 
Let be a positive Borel measure on M, with compact support. Its Cauchy transform 



F{z)^l- [ 

JR 



dfi{t) 



provides an analytic function on C \ R which is also regular at infinity, and has the 
normalizing value 1 there. The power expansion, for large values of \z\, yields the 
generating moment series of the measure /i: 



FU) = l + ^ + ^ + M^ + .... 



z'^ z^ 



On the other hand. 



lisiF{z) — —Imz J 



dn(t) 



whence 



lmF{z) Im^ < 0, z eC\R. 

Thus the main branch of the logarithm logF(^) exists in the upper half-plane and its 
imaginary part, equal to the argument of F{z), is bounded from below by — vr and from 
above by 0. According to Fatou's theorem, the non-tangential boundary limits 

f{t) = lim — Im logF(t + ie), 

exist and produce a measurable function with values in the interval [0, 1]. According to 
Riesz-Herglotz formula for the upper-half plane, we obtain: 



JR t — Z 



Or equivalently, 

F{z) = exp ' •'^ ' 



t 

One step further, let us consider the Lebesgue space L'^{n) and the bounded sclf- 
adjoint operator A = Mt of multiplication by the real variable. The vector ,^ = 1 
corresponding to the constant function 1 is 74-cyclic, and according to the spectral 
theorem: 



I 

jR 



t - z 

As a matter of fact an arbitrary function F which is analytic on the Riemann 
sphere minus a compact real segment, and which maps the upper/lower half-plane into 
the opposite half-plane has one of the above forms. These functions are known in 
rational approximation theory as Markov functions. 
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In short, putting together the above comments we can state the following result: 
the canonical representations: 



Fiz) 



exp( 



mdt. 

t-z' 



I -{{A 



t - z 

establish constructive equivalences between the following classes: 

a) Markov's functions F(z); 

b) Positive Borel measures /i of compact support on M; 

c) Functions / G L'^^^^{K) of compact support, < / < 1; 

d) Pairs (A, of bounded self-adjoint operators with a cyclic vector ^. 
The extremal solutions correspond, in each case exactly, to: 

a) Rational Markov functions F; 

b) Finitely many point masses /x; 

c) Characteristic functions / of finitely many intervals; 

d) Pairs (A, acting on a finite dimensional Hilbert space. 

For a complete proof see for instance Chapter VIII of |127] and the references cited 
there. The above dictionary is remarkable in many ways. Each of its terms has intrinsic 
values. They were long ago recognized in moment problems, rational approximation 
theory or perturbation theory of self-adjoint operators. 

For instance, when studying the change of the spectrum under a rank-one 
perturbation B = A — ^{■,^) one encounters the perturbation determinant: 

The above exponential representation leads to the phase-shift function fA,Bif) = f(t)'- 

fA,B{t)dt 



exp 



t - z 



The phase shift of, in general, a trace-class perturbation of a self-adjoint operator has 
certain invariance properties; it reflects by fine qualitative properties the nature of 
change in the spectrum. The theory of perturbation determinants and of the phase 
shift is nowadays well developed, mainly for its applications to quantum physics, see 

The reader will recognize above an analytic continuation in the complex plane of 
the real exponential transform 

f{t)dt- 



F{x) = Ef{x) = exp 



\t — x\ 



assuming for instance that x < M and the function / is supported by [M, oo). 

To give the simplest, yet essential, example, we consider a positive number r and 
the various representations of the function: 



Fiz) = 1 + 



z + r 



rdd^it) 



t 



exp 



dt 



t 



det[(— r — z){—zY 



In this case the underlying Hilbert space has dimension one and the two self-adjoint 
operators are A = Q and A — ^{■,^) = —r. 
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5.3.2. The reconstruction algorithm in one real variable Returning to our main theme, 
and as a direct continuation of the previous section, we are interested in the exact 
reconstruction of the original / : M — > [0, 1] from a finite set of its moments, or 
equivalently, from a Taylor polynomial of Ef at infinity. The algorithm described in 
this section is the diagonal Pade approximation of the exponential transform of the 
moment sequence. Its convergence, even beyond the real axis, is assured by a famous 
result discovered by A. A. Markov. 

Let oo, oi, . . . , a2n be a sequence of real numbers with the property that its 
exponential transform: 



exp 



C^O _|_ C^l _|_ 0,2n 



bo h 
1 + - + ^ + ..., 

Z Z'^ 



produces a non-ncgativc Hankcl matrix (bk+i)^i^Q. 

According to Markov's Theorem, there exists at least one bounded self-adjoint 
operator A G L{H), with a cyclic vector ^, such that: 



exp 



1 / Oq 0,1 0,2r. 



L\Z Z^ ^2n+l 



Let k <n and be the Hilbert subspace spanned by the vectors ^, A^, . . . , A^~^^. 
Suppose that dimi^^ = k, which is equivalent to saying that det(6i+j)fj'J^Q 7^ 0. Let tt^ 
be the orthogonal projection of H onto and let Ak — TTkATTk- Then 

whenever 0<i,j<A; — 1. In other terms, for large values of \z\: 

{{A - z)-'^, = {{Ak - z)-% + 

By construction, the vector ^ remains cyclic for the matrix e L(Hk). Let 
qk{z) be the minimal polynomial of A^, that is the monic polynomial of degree k which 
annihilates A^. In particular, 

qkiz){iAk- z)-^^,0 = {{qk{z) - qk{Ak)){Ak - z)-'i,0 = Pk-i{z) 

is a polynomial of degree k — 1. 
The two observations yield: 

q,(z){(A-z)-'^,0 = qk{z){{A,-z)-'^,0 + O{4^)=Pk-i{z) + O{- ^ ^ 



The resulting rational function Rk{z) — ^''g'^^^'' is characterized by the property: 



it is known as the Pade approximation of order {k — l,k), of the given series. 
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A basic observation is now in order: since ^'i; ■ ■ ■ ; ^2fe+i is the power moment 
sequence of a positive measure, qk is the associated orthogonal polynomial of degree k 
and pk is a second order orthogonal polynomial of degree A; — 1. In particular their roots 
are simple and interlaced. We prove only the first assertion, the second one being of 
a similar nature. Indeed, let ^ be the spectral measure of A localized at the vector ^. 
Then, for j < k, 



I 



t^qk(t)dt = {A^^,qk(A)0 = {Ai^,qk(Ak)0 = 0. 



Assume now that we are in the extremal case det(6i+j)"j-^o — that n is the 
smallest integer with this property, that is det(6j+j)"jiQ ^ 0. Since 

b,+j^{A%A^O, 

this means that the vectors ^, A$^, . . . , ^4"^ are linearly dependent. Or equivalently that 
Hn = H and consequently An = A. 

According to the dictionary established above, this is another proof that the 
extremal case of the truncated moment 1-problem with data ao, Oi, . . . , admits a 
single solution. The unique function / : R — > [0, 1] with this string of moments will 
then satisfy: 



exp 



n 



where the spectrum of the matrix A is {ai,...,an}, that of the perturbed matrix 
B — A — ^{-j^) is hi, ... ,bn and rj are positive numbers. Again, one can easily prove 
that 6i < Oi < &2 < 02 < . . . < &n < On- By the last example considered, we infer: 

n 

f ~ Xlbi,ai], 
1=1 

or equivalently 



. Pk-i + qk 
1 — sign 



Qk 

The above computations can therefore be put into a (robust) reconstruction 
algorithm of all extremal functions /. The Hilbert space method outlined above has 
other benefits, too. We illustrate them with a proof of another celebrated result due to 
A. A. Markov, and related to the convergence of the mentioned algorithm, in the case 
of non-extremal functions. 

Let n be a positive measure, compactly supported on the real line and let F{z) = 
Jjg(t — z)~^dfi{t) be its Cauchy transform. Then the diagonal Fade approximation 
Rn{z) = Pn-i{z) / qn{z) convcrgcs to F{z) uniformly on compact subsets o/C \ M. 
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This is the basic argument proving the statement: let A be the multiplication 
operator with the real variable on the Lebesgue space H = L'^{fi) and let ^ = 1 be its 
cyclic vector. The subspace generated by ^,A^, ...,^4""^^ will be denoted as before by 
Hn and the corresponding compression of Ahj An = TTnAiTn. 

If there exists an integer n such that H = Hn, then the discussion preceding the 
theorem shows that F = Rn and we have nothing else to prove. Assume the contrary, 
that is the measure n is not finite atomic. 

Let p{t) be a polynomial function, regarded as an element of H. Then 

{A - An)p{t) = tp{t) - (7r„A7r„)p(t) = tp{t) - tp{t) = 

provided that deg(p) < n. Since \\An\\ < \\A\\ for all n, and by Weierstrass Theorem, 
the polynomials are dense in H, we deduce: 

lim \\{A-An)h\\ =0, h e H. 

n— >oo 

Fix a point a G C \ M and a vector h E H. Then 

lim ||p-a)-i-(A„-a)->|| = lim \\{An - a)-\A - An){A - a)-^h\\ < 

n— >oo n— >oo 

1 



lim jT^MA - An){A - a)-^h\\ = 0. 



|Im a| 

A repeated use of the same argument shows that, for every /c > 0, 

lim IIP - a)-'' - {An - a)-'']h\\ = 0. 

n— >oo 

Choose a radius r < |Im a\ < \\{An — a)~^||~^, so that the Neumann series 

oo 

{An - z)-^ = {An-a-{z- a))-^ = '^{z - af{An - a) 



~k~l 
k=0 



converges uniformly and absolutely, in n and z, in the disk \z — a\ < r. Consequently, 
for a fixed vector h E H, 

lim \\{An - z)-^h -{A- z)-^h\\ = 0, 

n— >oo 

uniformly in z, \z — a\ < r. In particular, 

lim Rn{z) = {{An- z)-% = lim ((^n - z)-'^, = 

n— >oo n— >oo 

{{A-z)-XO = F{z), 

uniformly in z, \z — a\ < r. 

Details and a generalization of the above operator theory approach to Markov 
theorem can be found in |130] . 
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5.4- The exponential transform in two dimensions 



We return now to two real dimensions, and establish an analog of the matrix model for 
Markov's moment problem. Fortunately this is possible due to the import of some key 
results in the theory of semi-normal operators. We expose first the analog of Markov's 
exponential transform, and second, we will make a digression into semi-normal operator 
theory, with the aim at realizing the exponential transform in terms of (infinite) matrices, 
and ultimately of reconstructing planar shapes from their moments. 

The case of two real variables is special, partly due to the existence of a complex 
variable in . Let g : C — > [0, 1] be a measurable function and let dA{() stand for the 
Lebesgue area measure. The exponential transform of g, is by definition the transform: 

1 j-giOdAiC). 

TX 



Eg{z) = exp(-- 



This expression invites to consider a polarization in z: 



z E C \ supp((7). 



(180) 



Eg{z,w) = exp( 



9iC)dA{C) . 
^{C-z)iC-w)' 



z,w E C \ supp(5'). 



(6.1) 



The resulting function Eg{z, w) is analytic in z and antianalytic in w, outside the support 
of the function g. Note that the integral converges for every pair {z, w) G except the 
diagonal z = w. Moreover, assuming by convention exp(— oo) = 0, a simple application 
of Fatou's Theorem reveals that the function Eg{z,w) extends to the whole and it 
is separately continuous there. Details about these and other similar computations are 
contained in |127j . 

As before, the exponential transform contains, in its power expansion at infinity, 
the moments 



amn = amnig)= / z"" z"" g{z)dA{z) , m,n>0. 

According to Riesz Theorem these data determine g. We will denote the resulting series 
by: 



exp 



-1 



TT 



E 

m,n=0 



E 

m,n=0 



;i8i) 



The exponential transform of a uniformly distributed mass on a disk is simple, and 
in some sense special, this being the building block for more complicated domains. A 
direct elementary computation leads to the following formulas for the unit disk D, cf. 



Ej:)(z,w) 



zw 
z 



w 
w 



z, w G D 



z eB, w G D 



w eB, z eB 
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Remark that -Ed(^) = E-^^z, z) is a rational function and its value for |z| > 1 is 1 — 
The coefficients fo^n of the exponential transform are in this case particularly simple: 
6oo = 1 and all other values are zero. 

Once more, an additional structure of the exponential transform in two variables 
comes from operator theory. More specifically, for every measurable function g : 
C — ^ [0,1] of compact support there exists a unique irreducible, linear bounded 
operator T G L{H) acting on a Hilbert space if, with rank-one self-commutator 
[T*, T] = ^ ® ^ = ^(■, 0, which factors Eg as follows: 

Eg{z,w) = 1 - {{T* -w)-^i,{T* -zY^i), z,we supp((?)^ (182) 

As a matter of fact, with a proper extension of the definition of localized resolvent 
(T* — the above formula makes sense on the whole C^. The function g is called 

the principal function of the operator T. The next section will contain a brief incursion 
into this territory of operator theory. 

Let (yf : C — s> [0, 1] be a measurable function and let Eg{z,w) be its polarized 
exponential transform. We retain from the above discussion the fact that the kernel: 

1 — Eg{z,w), z,wEC, 

is positive definite. Therefore the distribution Hg{z,w) = —-^-^Eg{z,w) has compact 
support and it is positive definite, in the sense: 

/ Hg{z,w)(l){z)^{w)dA{z)dA{w) > 0, </> e C°°(C). 

If g is the characteristic function of a bounded domain Q G C, then it is elementary to 
see that the distribution Hq{z,w) = Hg{z,w) is given on f2 x by a smooth, jointly 
integrable function which is analytic in z G and antianalytic in w E fl, see |12Uj . 
In particular, this gives the useful representation: 

En{z, ti; = 1 - — / — f — — — — , z,w en , 

TT^ Jq2 {u- z){v -w) 

where the kernel Hq is positive definite in x fi. 

The example of the disk considered in this section suggests that the exterior 
exponential transform of a bounded domain E^{z,w) may extend analytically in each 
variable inside fl. This is true whenever dfl is real analytic smooth. In this case there 
exists an analytic function S defined in a neighborhood of dQ, with the property: 

S{z) = z, ze on. 



The anticonformal local refiection with respect to dQ is then the map z ^ S{z); for this 
reason S{z) is called the Schwarz function of the real analytic curve dQ, introduced 
earlier in this text. Let a; be a relatively compact subdomain of Q, with smooth 
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boundary, too, and such that the Schwarz function S{z) is defined on a neighborhood 
of f2 \ u;. A formal use of Stokes' Theorem yields: 



But the latter integral is analytic/ antianalytic for z,w & uJ^. A little more work with 
the above Cauchy integrals leads to the following remarkable formula for the analytic 
extension of Eq{z,w) from z,w & if to z,w & u'^: 



The study outlined above of the analytic continuation phenomenon of the exponential 
transform Eq{z, w) led to a proof of a priori regularity of boundaries of domains which 
admit analytic continuation of their Cauchy transform. The most general result of this 
type was obtained by different means by Sakai. We simply state the result, giving in this 
way a little more insight into the proof of the regularity of the boundaries of quadrature 
domains. 

Let Q be a bounded planar domain with the property that its Cauchy transform 



extends analytically across dQ. Then the boundary dQ is real analytic. 

Moreover, Sakai has classified the possible singular points of the boundary of such 
a domain. For instance angles not equal to or vr cannot occur on the boundary. 

5.5. Semi-normal operators 

A normal operator is modelled via the spectral theorem as multiplication by the complex 
variable on a vector valued Lebesgue L^-space. The interplay between measure theory 
and the structure of normal operators is well known and widely used in applications. 
One step further, there are by now well understood functional models, and a complete 
classification for classes of close to normal operators. We record below a few aspects 
of the theory of semi-normal operators with trace class self-commutators. They will 
be serve as Hilbert space counterparts for the study of moving boundaries in two 
dimensions. The reader is advised to consult the monographs \127\ 1131] for full details. 

Let if be a separable, complex Hilbert space and let T G C{H) be a linear 
bounded operator. We assume that the self-commutator [T*,T] = T*T — TT* is trace- 
class, and call T semi-normal. If [T*,T] > 0, then T is called hypo-normal. For 
a pair of polynomials p{z,z),q{z,z) one can choose (at random) an ordering in the 






{z — S{w)){S{z) — w)Hq{z, w), z,w E Q\u!. 
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functional calculus p{T,T*),q{T,T*), for instance putting all adjoins to the left of all 
other monomials. The functional 

(p,g)^trace[p(T,T*),g(T,r*)] 

is then well-defined, independent of the ordering in the functional calculus, and possesses 
the algebraic identities of the Jacobian A direct (algebraic) reasoning will imply 

the existence of a distribution ut G ^^'(C) satisfying 

tTace[p{T,T*),qiT,T*)]=UT 

see |132] . The distribution ut exists in any number of variables (that is for tuples of 
self-adjoint operators subject to a trace class multi-commutator condition) and it is 
known as the Helton-Howe functional. 

Dimension two is special because of a theorem of J. D. Pincus which asserts that 
'^T = ^grdA, that is ut is given by an integrable function function qt, called the 
principal function of the operator T, see [133[ 1134] . 

The analogy between the principal function and the phase shift (the density of the 
measure appearing in Markov's moment problem in one variable) is worth mentioning 
in more detail. More precisely, if B = A — K is a. trace-class, self-adjoint perturbation 
of a bounded self-adjoint operator A G L{H), then for every polynomial p{z), Krein's 
trace formula holds: 

tT[piB)-piA)]= [ p'it)fAAt)dt, 

where /a,s is the corresponding phase-shift function, |128j . It is exactly this link between 
Hilbert space operations and functional expressions which bring the two scenarios very 
close. Taking one step further, exactly as in the one variable case, the moments of the 
principal function can be interpreted in terms of the Hilbert space realization, as follows: 

mk j z'^-^z^-^gT{z)dA = 

trace [T*^T'"], k,m>l. 

In general, the principal function can be regarded as a generalized Fredholm index of 
T, that is, when the left hand side below is well defined, we have 

ind{T-X) = -gT{X). 

Moreover enjoys the functoriality properties of the index, and it is obviously invariant 
under trace class perturbations of T. Moreover, in the case of a fully non-normal 
operator T, 

supp^fT = cr(T), 

and various parts of the spectrum cr{T) can be interpreted in terms of the behavior of 
gr, see for details |127] . 



d{p, q) 
d{z, z) 
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{Tfm=tf{t)+^[u{t)f{t) 



To give a simple, yet non-trivial, example we proceed as follows. Let be a planar 
domain bounded by a smooth Jordan curve F. Let H'^{T) be the closure of complex 
polynomials in the space L^{T, ds), where ds stands for the arc length measure along F 
(the so-called Hardy space attached to F). The elements of -ff^(F) extend analytically 
to Q. The multiplication operator by the complex variable, T^f = zf, f G if^(F), 
is obviously linear and bounded. The regularity assumption on F implies that the 
commutator [T^, T*] is trace class. Moreover, the associated principal function is the 
characteristic function of Q, so that the trace formula above becomes: 

trace[p(T„T;),g(T„T;)] = - / ^|?^dA, p,qeC[z,z]. 

See for details [T271[T8T] . 

A second, more interesting (generic example this time) can be constructed as 
follows. Let u{t),v{t) be real valued, bounded continuous functions on the interval [0, 1]. 
Consider the singular integral operator, acting on the Lebesgue space L'^{[0,l],dt) by 
the formula: 

1 v{t)v{s)f{s)ds 
■n Jq s-t 

Then it is easy to see that the self-commutator [T*,T] is rank one. The principal 
function qt will be in this case the characteristic function of the closure of the domain 
G given by the constraints 

G = {{x,y) G M^; \y - u{x)\ < ^;(a;)^ x G [0,1]}. 

Based on a refinement of this example, in general every hyponormal operator with trace 
class self-commutator can be represented by such a singular integral model, with matrix 
valued functions u, v, acting on a direct integral of Hilbert spaces over [0, 1] ; in which case 
the principal function relates directly to Krein's phase shift, by the following remarkable 
formula due to Pincus |133j : 

5't(2^; 1/) fu{x)—v{x)*v(x),u{x)+v{x)*v{x){y)- 

The case of rank-one self-commutators is singled out in the following key 
classification result: 

There exists a bijective correspondence T ^ qt between irreducible hyponormal 
operators T , with rank-one self-commutator, and bounded measurable functions with 
compact support in the complex plane. 

An invariant formula, relating the moments of the principal function g to the 
Hilbert space operator T, [T*,T] = ^(■,^),, satisfying qt = g,a.e. is furnished by the 
determinantal formula: 

exp(-- / -/^^^^^^) = det[(T*-«;)-i(T-^)(r*-Iz7)(T-.)-i] = 
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This formula explains the positivity property of the exponential transform, alluded to 
in the previous section. 

The bijective correspondence between classes g G L^j^p(C), < 5^ < 1 and 
irreducible operators T with rank-one self-commutator was exploited in |135l I136j for 
solving the L-problem of moments in two variables. The theory of the principal function 
has inspired and played a basic role in the foundations of modern non-commutative 
geometry (specifically the cyclic cohomology of operator algebras) and non- commutative 
probability. 

We have to stress the fact that the above bijective correspondence between "shade 
functions" Qt and irreducible hyponormal operators T with rank-one self-commutator 
can in principle transfer any dynamic g{t) into a Hilbert space operator dynamic T{t). 
However, the details of the evolution law of T{t) even in the case of elliptic growth 
are not trivial, nor make the integration simpler. We will see some relevant low degree 
examples in the next section. 

5.5.1. Applications: Laplacian growth To give a single abstract illustration, consider a 
growing family of bounded planar domains D{t) with smooth boundary: 

D{t) C D{s), whenever t < s. 

The evolution of the exponential transforms 

-1 f dAiC) 



En(t)iz,w) = exp 



is governed by the differential equation (in the standard vector calculus notation) 

-EDit){z,w) = —ED(t){z,w) - . 

dt vr Jg^^^t) (( - z){C- w) 

Any evolution law at the level of the pair (T(t),^(t)) will have the form 

^^Enit){z,w) = {{T*{t) -wr')T*{t){T*{t) -w)-'mdT%t) -z)-'m)- 
{{T*{t)-w)-'e{t),{T*{t)-z)-'m)+ 

((T*(t) - wr'm, (T*(t) - zr'T*{t)iT*it) - z^m)- 
{{T*{t)-wr'm,{T*{t)-z)-')^'{t)). 

A series of simplification in the case of elliptic growth are immediate: for instance 
is proportional to the area of D{t), whence we can choose the vector of the form 

Second, the higher harmonic moments are preserved by the evolution, whence the 
Cauchy transform/resolvent 

^^m,iT*it)-z)-'m)- ^ '^'^^^^ - ' 



dt dt Jd^i) C ~ z z 
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gives full information about the first row and first column in the matrix representation of 
T*{t) in the basis obtained by orthonormalizing the sequence T*(t)^(t), T*^^(t), .... 
The reader can consult the article [122] for more details about computations related to 
the above ones. 



5.6. Linear analysis of quadrature domains 

If we would infer from the one-variable picture a good class of extremal domains for 
Markov's L-problem in two variables we would choose the disjoint unions of disks, as 
immediate analogs of disjoint unions of intervals. In reality, the nature of the complex 
plane is much more comphcated, but again, fortunately for our survey, the class of 
quadrature domains plays the role of extremal solutions in two real dimensions. 

Recall from our previous sections that a bounded domain Q of the complex plane 
is called a quadrature domain (always henceforth for analytic functions) if there exists a 
finite set of points Oi, 02, . . . , G fl, and real weights Ci, C2, . . . , Q, with the property: 

[ f{z)dA{z) = ci/(ai) + C2/(a2) + . . . + Q/(a,), / G AL\Q) 
Jn 

where the latter denotes the space of all integrable analytic functions in Q. In case some 
of the above points coincide, a derivative of / can correspondingly be evaluated. 
Let Q he a bounded planar domain with moments 

aran = amni^)= / z"'z''dA{z), m, n > 0. 
Jn 

The exponential transform produces the sequence of numbers bmn = bmn{^), m,n > 0. 
Let T denote the irreducible hyponormal operator with rank-one self-commutator 



[T*,r] = ^{■,0- In virtue of the factorization (182) 



fe,nn = (T*™e,r*'^e), m,n>0. 

Hence the matrix {bmn)mn=o turns out to be non-negative definite. The following result 
identifies a part of the extremal solutions of the L-problem of moments as the class of 
quadrature domains: 

A bounded planar domain Q is a quadrature domain if and only if there exists a 
positive integer d > 1 with the property det(6mn(^))m,n=o — 0- 

For a proof see |135] . The vanishing condition in the statement is equivalent to the 
fact that the span Hd of the vectors ^, T*^, T*^^, ... is finite dimensional (in the Hilbert 
space where the associated hyponormal operator T acts). Thus, if f2 is a quadrature 
domain with corresponding hyponormal operator T, and is the compression of T to 
the (i-dimensional subspace H^, then: 

Eniz, w) = l- ((T; - w)-'^, {T: - z)-^i), z,we 



In particular this proves that the exponential transform of a quadrature domain is a 
rational function. As a matter of fact a more precise statement can easily be deduced: 
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Let n be the quadrature domain defined above. Then 

^ ' P(z)P(w) 



This result offers an efficient characterization of quadrature domains in terms of a 
finite set of their moments (see the reconstruction section below) and it opens a natural 
correspondence between quadrature domains and certain classes of finite rank matrices. 
We only describe a few results in this direction. For more details see |120l 11221 1135j . 

In the conditions of the above result, let f2 be a quadrature domain with associated 
hyponormal operator T; let Hq = Vfc>o T*''^ and let p denote the orthogonal projection 
of the Hilbert space H (where T acts) onto Hq. Denote Cq = pTp (the compression of 
T to the d-dimensional space Hq) and Dq = [T*,T]. Then the operator T has a two 
block-diagonal structure: 

\ 



/ Co 











Di 


























Ds 





V 



where the entries are a\\ d x d matrices, recurrently defined by the system of equations: 

[Ck*, Ck] + Dk+i* Dk+i = DkDk* 
Cfc+i*-Dfc+i = Dk+iCk* , k > 0. 

Note that > for all k. This decomposition has an array of consequences: 

(i) The spectrum of Cq coincides with the quadrature nodes of Q; 

(ii) n = {z; II (Cq - ^)"^^|| > 1} (up to a finite set); 

(iii) The quadrature identity becomes 



/(^)dA(^) = vr(/(Co)e,e), 



for / analytic in a neighborhood of Q; 
(iv) The Schwarz function of Q is 



S{z) 



where z E Q. 



To give the simplest and most important example, let $7 = D be the unit disk (which 
is a quadrature domain of order one) . Then the associated operator is the unilateral 
shift T = Tz acting on the Hardy space H'^{dT)). Denoting by z"' the orthonormal basis 
of this space we have Tz"' = z"'~^^, n > 0, and [T*,T] = 1) is the projection onto 
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the first coordinate 1 = z^. The space Hq is one dimensional and Cq = 0. This will 
propagate to = and D^. = 1 for all k. Thus the matricial decomposition of T 
becomes the familiar realization of the shift as an infinite Jordan block. 

In view of the linear algebra realization outlined in the preceding section we obtain 
more information about the defining equation of the quadrature domain. For instance: 

Q{z,w) 



P{z)P{w) 



which yields 



i-l 



g(z,^) = |p(^)p-5^|g,(2)p, 

fc=0 

where is a polynomial of degree k in z, see |122] . 

Thus the exponential transform of a quadrature domain contains explicitly the 
irreducible polynomial Q which defines the boundary and the polynomial P which 
vanishes at the quadrature nodes. By putting together all these remarks we obtain 
a strikingly similar picture to that of a single variable. More specifically, if f2 is a 
quadrature domain with d nodes, as given above, and associated hyponormal operator 
T, then: 

Eniz,w) = -^^^ = 1 - {iT:-w)-'^,iT:--z)-'0 = 
P[z)P[w) 

d _ 



EHn{ai,aj) — _ ^ _ , z.weil" 
n ■ — 7' n ■ — in 



TT^ — ' " ai — z an — w 

i,j=i -I 



In particular we infer, assuming that all nodes are simple: 

~^ o,/ ^-l^7^ = CiCjHn{a,,aj), l<t,j <d. 
P'{ai)P\aj) 

For details see [T35|[T22]. 

The interplay between these additive, multiplicative and Hilbert space decomposi- 
tions of the exponential transform gives an exact reconstruction algorithm of a quadra- 
ture domain from its moments. The next section will be devoted to this algorithm. 

Before ending the present section we consider an illustration of the above formulas. 
Let VL = uf^iD{ai,ri) be a union of d pairwise disjoint disks. This is a quadrature 
domain with data: 

P{z) = {z - ai) . . . {z - ad), 

Q{z,w) = [{z - ai)(w -ai) -rf]... [{z - ad){w -Od) - r'^dl- 

The associated matrix is also computable, involving a sequence of square roots of 
matrices, but we do not need here its precise form. Whence the exponential transform 
is, for large values of \z\, \ w\: 



rf 1 1 , Q{ai,aj) Vi r 



- - fli) ijr^i P'{ai)P'{aj) ai~zaj-w' 
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The essential positive definiteness of the exponential transform of an arbitrary 
domain can be deduced, via an approximation argument, from the positivity of the 
matrix {—Q{ai,aj))fj^-j^, where Q is the defining equation of a disjoint union of 
disks. We note that (— a7))i!i=i > is only a necessary condition for the disks 
D{ai,ri), 1 < i < d, to he disjoint. Exact computations for d = 2 immediately show 
that this matrix can remain positive definite even the two disks overlap a little. However, 
if two disks overlap, then, by adding an external disk, even far away, this prevents the 
new 3x3 matrix to be positive definite. 

We end this section with two examples, covering the totality of quadrature domains 
of order two. 

Quadrature domains with a double node. Let z = w"^ + bw he the conformal 
mapping of the disk \w\ < 1, where b > 2. Then z describes a quadrature domain Q of 
order 2, whose boundary has the equation: 

Q{z,z) = \z\^-{2 + b^)\z\^ -b^z-b''z + l-b^ = 0. 

The Schwarz function of Q has a double pole at z = 0, whence the associated 
2 X 2-matrix Cq is nilpotent. Moreover, we know that: 

\z\'UC*-z)-'a' = \z\'-P{z,z). 

Therefore 

||(C* + z)^\\^ = (2 + b^)\z\^ + b^z + b'^z + 6^-1, 

or equivalently: \\if = 2 + 6^ (C*e,0 = and \\Clif = b'^-l. 
Consequently the linear data of the quadrature domain VL are: 

/ ^^-1 \ / —Jl 

r* — \ (fe2-2)l/2 1 C — \ (b2-l)l/2 

Quadrature domains with two distinct nodes. Assume that the nodes are 
fixed at ±1. Hence P{z) = z"^ — 1. The defining equation of the quadrature domain Q 
of order two with these nodes is: 

Q{z,z) = i\z + l\^-r^){\z~l\^-r^)-c, 

where r is a positive constant and c > is chosen so that either is a union of two 
disjoint open disks (in which case c = 0), or Q{0, 0) = 0, see [109] . A short computation 
yields: 

Q{z, z) = z^z^ — 2rzz — z'^ — z'^ + cn^r), 

where 

r < 1 
r > 1. 



air) 




One step further, we can identify the linear data from the identity: 

\P{zni-\\{C;-z)-'a') = Q{z,z). (183) 
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Consequently, 

^ \ u / \ y 

\f2r 

This simple computation illustrates the fact that, although the process is affine in 
r, the linear data of the growing domains have discontinuous derivatives at the exact 
moment when the connectivity changes. 

5. 7. Signed measures, instability, uniqueness 

Contrary to the uniqueness of a quadrature domain for subharmonic functions with a 
prescribed quadrature measure, quadrature domains for harmonic or analytic functions 
are not determined by the quadrature nodes and weights. This is an intriguing global 
phenomenon which has haunted mathematicians for many decades, we briefly record 
below some significant discoveries in this direction. 

Consider quadrature domains for harmonic test functions and real- valued measures 



(174). As to the relationship between the geometry of and the location of supp fi there 
are then drastic differences between the cases of having all Cj > respectively having 
no restrictions on the signs of Cj. This is clearly demonstrated in the following theorem 
due to M. Sakai |137] , |138] . The second part of the theorem is discussed (and proved) 
in some other forms also in |lU8j . |139j . |1U6] . |1U7] . |14U] . for example: 



Let r and R be positive numbers, R > 2r. Consider measures (i of the form (174) 
with Cj real and related to r and R by 

supp/i C 5(0,r), (184) 

n 

^c, = rf. (185) 

j=i 

(i) If fi > 0, then any quadrature domain Q for harmonic functions for fi is also a 
quadrature domain for subharmonic functions. Hence the previous result applies, 
and in addition 

B{0,R-r) cnc B{0,R + r). 

(a) With fi not necessarily > 0, and with no restrictions on ^2^=1 kjl ond n, 
any bounded domain containing B{0,r) and having area ttR"^ can be uniformly 
approximated by quadrature domains for harmonic functions for measures /i 



satisfying (184), (185) 



With /i a signed measure of the form (174) we still have YTj=i = 1^1 5 but ^"=1 



may be much larger. In view of the theorem, the ratio 

^ E?=i Cj ^ ^ dfi 
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(0 < p < 1) might give an indication of liow strong is the coupling between the geometry 
of supp fi and the geometry of Q. 

As mentioned, a quadrature domain for harmonic functions is not always uniquely 
determined by its measure /i. Still there is uniqueness at the infinitesimal level: if 

n „ 

y^cM^j) = hpdA (186) 

j=i 

and (for example) the aj are kept fixed, then one can always increase the Cj (indefinitely) 
and get a unique evolution of Q (Hele-Shaw evolution). If dfl has no singularities then 
one can also decrease the cj slightly and have a unique evolution (backward Hele-Shaw, 
which is ill-posed). Thus it makes sense to write 

n = i7(Ci, . . . , Cn) 

for Cj in some interval around the original values. Note however that decreasing the Cj 
makes the ratio p decrease, indicating a loss of control or stability. 

In the simply connected case, Q will be the image of the unit disc D under a rational 
conformal map / = /{ci,...,c„) : D ^{ci, . . . , c„). This rational function is simply the 
conformal pull-back of the meromorphic function {z,S{z)) on the Schottky double of 
Q to the Schottky double of D, the latter being identified with the Riemann sphere. 
It follows that the poles of / are the mirror points (with respect to the unit circle) 
of the points /"^(a^). When the Cj increase then the \f~^{aj)\ decrease (this follows 
by an application of Schwarz' lemma to /i~i.gerc ° /originaicj), hence the poles of / move 
away from the unit circle. Conversely, the poles of / approach the unit circle as the Cj 
decrease, also indicating a loss of stability. 

For decreasing Cj the evolution fi(ci,...,c„) always breaks down by singularity 
development of dQ or dQ reaching some of the points aj (see e.g. |1U3] . |141] ) before fl 
is empty, except in the case that Q{ci, . . . , c„) is a quadrature domain for subharmonic 
functions. In the latter case the Cj (necessarily positive) can be decreased down to zero, 
and fl will be empty in the limit ci = . . . = c„ = 0. However, it may happen that 
fi(ci, . . . , c„) breaks up into components under the evolution. 

Assume now that Q is simply connected. Then the analytic and harmonic functions 



are equivalent as test classes for (186). In the limit case that all the points Uj coincide. 



say ai = . . . = a„ = 0, then (186) corresponds to 



n 



^cy^--i)(0)= / vdA (187) 

for if analytic. The cj (allowed to be complex) now have a slightly different meaning 
than before. In fact, they are essentially the analytic moments of fl: 



Cj = ^. ^ j ^dA {j = l,...,n) 



The higher order moments vanish, and the conformal map / = f{ci,...,cn) '■ D 
fi(ci,...,c„) (normalized by /(O) = 0, /'(O) > 0) is a polynomial of degree n. A 
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precise form of the local bijectivity of the map (ci, . . . , c„) ^{ci, . . . , c„) has been 
established by O. Kouznetsova and V. Tkachev |142j . |143j . who proved an explicit 
formula for the (nonzero) Jacobi determinant of the map from the coefficients of / to 
the moments (ci, . . . ,c„)). This formula was conjectured (and proved in some special 
cases) by C. Ullemar |144] . 



On the global level, it does not seem to be known whether (187), or (186), with a 
given left member, can hold for two different simply connected domains and all analytic 
If. 

Leaving the realm of quadrature domains, an explicit example of two different 
simply connected domains having the same analytic moments has been given by M. Sakai 
[145] . The idea of the example is that a disc and a concentric annulus of the same area 
have equal moments. If the disc and annulus are not concentric, then the union of them 
(if disjoint) will have the same moments as the domain obtained by interchanging their 
roles. Arranging everything carefully, with removing and adding some common parts, 
two different Jordan domains having equal analytic moments can be obtained. Similar 
examples were known earlier by A. Celmins |146j . and probably even by P. S. Novikov. 
On the positive side, a classical theorem of Novikov |147j asserts that domains which 
are starshaped with respect to one and the same point are uniquely determined by their 
moments. See |148] for further discussions. 

Returning now to quadrature domains, there is definitely no uniqueness for 
harmonic and analytic test classes if multiply connected domains are allowed. If Q 
has connectivity m + 1 (m > 1), i.e., has m "holes", then there is generically an m- 
parameter family fi(ti, . . . , t^) of domains such that fl{0, . . . ,0) = Q and 



_d_ 



ipdA = (j = l,...,m) 

n{ti , . . . ,tm) 



for every if analytic in a neighborhood of the domains. These deformations are Hele- 
Shaw evolutions, driven not by Green functions but by "harmonic measures" , i.e., regular 
harmonic functions which take (different) constant boundary values on the components 

of dn. 

It follows that multiply connected quadrature domains for analytic functions for a 
given /i occur in continuous families. It even turns out |18], |149] that any two algebraic 
domains for the same n can be deformed into each other through families as above. 
Thus there is a kind of uniqueness at a higher level: given any /i there is at most one 
connected family of algebraic domains belonging to it. 

For harmonic quadrature domains there are no such continuous families (choosing 



(p{z) = \og\z — a\ in (186) with a E C \ ^7 in the holes stops them), but one 
can still construct examples with a discrete set of different domains for the same /x. 
It is for example possible to imitate the example with a disc and an annulus with 
quadrature domains for measures /i of the form (174), with > 3) and 



ci = . . . = c„ = c > suitably chosen. However, it seems very difficult to imitate 
the full Sakai construction, with "removing and adding some common parts", in the 
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context of quadrature domains. Therefore it is not at all easy to construct different 
simply connected quadrature domains for the same /i. 

We end this section with the simplest example of a continuous class of quadrature 
domains with the same quadrature data. 

Three points, non-simply connected quadrature domains and the non- 
uniqueness phenomenon. Quadrature domains (for analytic functions) with at most 
two nodes, as in the above examples, are uniquely determined by their quadrature data 
and are simply connected. For three nodes and more it is no longer so. The following 
example, taken from |lU9j . with three nodes and symmetry under rotations by 27r/3, 
illustrates the general situation quite well. More details on the present example are 
given in |1U9] . and similar examples with more nodes are studied in |121] . 

Let the quadrature nodes and weights be Oj = and Cj = vrr^ respectively 
(j = 1,2,3), where u = e^'^*/^ and where r > is a parameter. Considering first the 



strongest form of quadrature property, namely for subharmonic functions, as in (173) 



(174), the situation is in principle easy: f2 is for any given r > uniquely determined 



up to nuUsets and can be viewed as a swept out version of the quadrature measure 
H = Cj<5aj or as the union of the discs B{aj,r) with (possible) multiple coverings 

smashed out. 

For < r < ^ the above discs are disjoint, hence f2 = U^^iB{aj,r). For r larger 
than ^ but smaller than a certain critical value ro (which seems to be difficult to 
determine explicitly) Q is doubly connected with a hole containing the origin, while 
for r > ro the hole will be filled in so that f2 is a simply connected domain. The 
above quadrature domains (or open sets) are actually uniquely determined even within 



nuUsets, except in the case r = tq when both f2 and \ {0} satisfy (173). 

Consider next the general class of quadrature domains for analytic functions 
(algebraic domains). For < r < ^ only the disjoint discs qualify, as before. However, 
for any r > ^^ there is a whole one-parameter family of domains satisfying the 
quadrature identity for analytic ip. These are defined by the polynomials 

Q{z, z) = z^z^ - z^ -z^ - Sr^z^z'^ - (188) 

3r(r^ - 2rV + l)zz + T^{2r^ - 3rV + 1), 

where r > is a free parameter, independent of the quadrature data. When completed 
as to nullsets, the quadrature domains in question are more precisely 

Q{r,T) = intclos{z G C : Q{z,z) < 0}. 

The interpretation of the parameter r is that on each radius {z = tuj^^^ : t > 0}, 
j = 1,2,3, there is exactly one singular point of the algebraic curve Q{z,z) = 0, and 
r = \z\ for that point. This singular point is either a cusp on dQ or an isolated 
point of Q{z,z) = 0, a so-called special point. Special points are those points a G f2 for 
which the quadrature identity admits the (integrable) meromorphic function (p{z) = 
Equivalently, f2\{a} remains to be a quadrature domain for integrable analytic functions. 
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For ^ < r < the quadrature domains for analytic functions are exactly 
the domains fi(r, r) (with possible removal of special points) for r in an interval 
< T < T2(^), where ri(r), T2(r) satisfy < ri(r) < ^ < r2(r), and more precisely 
can be defined as the positive zeros of the polynomial 4r^ — 4r^r + 1. (see |109] for 
further explanations and proofs). The domains fl{r, r) are doubly connected with a hole 
containing the origin. When r increases the hole shrinks and both boundary components 
move towards the origin. For r = T2(r) there are three cusps on the outer boundary 
component which stop further shrinking of the hole, and for r = ri(r) there are three 
cusps on the inner boundary component which stop the expansion of the hole. 

For exactly one parameter value, r = Tg^hhi^), ^('"jt) is a quadrature domain 
for subharmonic functions (and so also for harmonic functions). This Tsubh('") can be 
determined implicitly by evaluating the quadrature identity for ip{z) = log|2;|, which 
gives the equation 

/ log|;z|dA(;z) = 0. 

For r = ri(r) = T2{r) = |, and as r increases, Ti(r) decreases and T2(r) 
increases. What happens when r = 2~6 is that for f2(r, T2(r)), i.e., for the domain with 
cusps on the outer component, the hole has shrunk to a point (the origin). Hence, for 
r = 2~s, n{r,T2{r)) is simply connected, while f2(r, r) for ri(r) < r < T2(r) remain 
doubly connected. 

For all ^ < r < 2~^, ri(r) < rsubh('") < T2{r) because a subharmonic quadrature 
domain cannot have the type of cusps which appear for r = ri(r), r2(r) (see |99], |100j ). 
It follows that the critical value r = tq, when f2(r, rsubh('")) becomes simply connected, 
is larger that 2~i. 

For r > 2~6 the quadrature domains for analytic functions are the domains fi(r, r) 
(with possible deletion of special points), with r in an interval Ti(r) < r < Ts{r). Here 
ri(r) is the same as before (i.e., corresponds to cusps on the inner boundary), while T^{r) 
is the value of r for which the hole at the origin degenerates to just the origin itself 
(which for r > 2^^ occurs before cusps have developed on the outer boundary). The 
origin then is a special point, and one concludes from (188) that r = r-^ir) is the smallest 
positive zero of the polynomial 2r^ — 3r^r + 1. For r = 2~5, r^ir) = r2(r) = 2~5. 

For 2~6 < r < To we have Ti(r) < Tsuhhi^) < r^^r), while for r > Tq, Tsubh('") = ''"3('")- 
Since Q{r,T3{r)) is simply connected and is a quadrature domain for analytic functions 
it is also a quadrature domain for harmonic functions. It follows that in the interval 
2~g < r < tq there are (for each r) two different quadrature domains for harmonic 
functions, namely ^{r,Tsuhh{r)) and fl{r,T3{r)) (doubly respectively simply connected). 

In summary, we have for each r > ^ a one-parameter family of algebraic domains 
f2(r, r), for exactly one parameter value (r = rsubh('")) this is a quadrature domain for 
subharmonic functions, and for each r in a certain interval (2^i < t < rg) there are two 
different quadrature domains for harmonic functions (f^(r, rsubh('")) and Q{r,T3{r))). 
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6. Other physical appUcations of the operator theory formulation 

The preceding chapters provide a review of the relationships between the theory of 
normal random matrices, where evolution is defined by increasing the size of the matrix 
(a discrete time), its continuum (or infinite size) limit - Laplacian growth - and the 
general theory of semi-normal operators whose spectrum approximates generic domains. 
The exposition refiects, to some extent, the parallel historical development of the two 
non-commutative generalizations of Laplacian growth (random matrix theory and semi- 
normal operator theory). It is quite natural, at this point, to investigate the direct 
relationships between these two theories. However, this is a task of a magnitude which 
would require a separate review at the very least. We will therefore contend ourselves 
with exposing only a few of these relations, via their applications to physical problems. 

The first application has to do with refined asymptotic expansions which 
characterize Laplacian growth in the critical case, before formation of a (2, 3) cusp. 
As we will see, to obtain this limit, one must take a "double-scahng limit" by fine- 
tuning two parameters of the random matrix ensemble. Alternatively, this procedure is 
equivalent to a special choice of Pade approximants in the operator theory approach. 

The second application described in this section is a very brief introduction of the 
notion of free, non- commutative random variables, and its relevance in open problems of 
strongly interacting quantum models, particularly in the 2D metal-insulator transition 
and the determination of ground state for 2D spin models. The review concludes with 
this cursory exposition. 

6.1. Cusps in Laplacian growth: Painleve equations 

In this section, we exploit the formalism built up to now, in order to address a problem 
of great significance both at the mathematical and physical levels: what happens when 
a planar domain evolving under Laplacian growth approaches a generic (2, 3) cusp? 
We have already seen that a classical solution does not exist, in that no singly- 
connected domain with uniform density would satisfy the conditions of the problem. 
However, since we now have alternative formulations of Laplacian growth via the 
balayage of the uniform measure, we may generalize the problem and ask whether there 
is any equilibrium measure, dropping the uniformity (and indeed, the two-dimensional 
support) of the classical solution. By analogy with the ID situation, we seek a solution 
in the sense of Saff and Totik, where the support and density of the equilibrium measure 
are given by the proper weighted limit of orthogonal polynomials. In order to obtain this 
limit, we must organize the evolution equations of the wavefunction such as to extract 
the correct scaling limit, for — > oo. 

6.1.1. Universality in the scaling region at critical points - a conjecture Detailed 
analysis of critical Hermitian ensembles indicates that the behavior of orthogonal 
polynomials in a specific region including the critical point (the scaling region), upon 
appropriate scaling of the degree n, is essentially independent of the bulk features of the 
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ensemble. This universality property (a common working hypothesis in the physics of 
critical phenomena) is expected to occur for critical NRM ensembles as well - and is 
indeed easy to verify in critical Gaussian models, 2|t2| = 1- Analytically, it means that 
by suitable scaling of the variables z,n: 

n oo, ^ 0, nh = to, to = tc~ fi^^^i z = Z(. + ffC,i 

where is the location of the critical point and t^ the critical area, the wave function 
'^n{z) will reveal a universal part (/>(z^, Q which depends exclusively on the local singular 
geometry ~ y'^ {p, q mutual primes) of the complex curve at the critical point. This 
conjecture is a subject of active research. Its main consequence is that in order to 
describe the scaling behavior for a certain choice of p, q, it is possible to replace a given 
ensemble with another which leads to the same type of critical point, though they may 
be very different at other length scales. 

6.1.2. Scaling at critical points of normal matrix ensembles In the remainder of the 
section we analyze the regularization of Laplacian Growth for a critical point of type 
p = 3, g = 2, by discretization of the conformal map as described in the previous 
paragraph. For simplicity, we start from the conformal map corresponding to the 
potential V{z) = t^z^, which is the simplest model leading to the specified type of 
cusp. It should be noted that the analysis will be identical for any monomial potential 
V{z) = tnZ"',n > 3; for every such map, n singular points of type p = 3,q = 2 will 
form simultaneously on the boundary. The critical boundary corresponding to n = 3 is 
shown in Figure [7} 

The scaling limit from the string equation We start from the Lax pair corresponding 
to the potential V{z) = t^z^: 

Li)n = Tnlpn+l + -^^Vn = Tn-llpn-l + M„+2^^n+2- (189) 



The string equation {'bl ) [L\ L\ = h translates into 



{'^1 + \u\n? - rl_i - \u\l^2)'>Pn + (r^Un+S " ?^n+2Mn+2)^n+3 

+ {rn~3Un - rn~lUn-l)iJn~3 = ^^n- (190) 

Identifying the coefficients gives 

{rl - |m|^+2 - - {rl^i - \u\l+i - \u\l) = h, (191) 



and 



3t3. (192) 



Equation (191) gives the quantum area formula 

- {\u\l+2 + = nh, (193) 



which together with the conservation law (192) leads to the discrete Painleve equation 

rl[l-9\ts\'{rl,,+rl^,)]=nh. (194) 
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In the continuum limit, the equation becomes 

- 18|t3|V = io- 
The critical (maximal) area is given by 



0, S6\h\''rl = 1. 



Choosing Tc = 1 gives 6|t3| ~ 1 and tc = \. It also follows that 



_ rn-2rn-i _ 3 
Un- 2 ' 2' 



Introduce the notations 



Nh = tc, nh = to = tc + n^V, rl = l- /^^"^(i/), z = - + h^X, 



where a = |- We get — tf-d^ and 

1 - n^^u - kh^^'uiu) 



n+k 



(195) 
(196) 

(197) 

(198) 

(199) 



where dot signifies derivative with respect to u. The scaling limit of the quantum area 
formula becomes 



(1 - h^^'u) 

giving at order h'^"' the Painleve 1 equation 

u-2u^ = 4i/. 
Rescaling u — > C2U, v — > c\v gives the standard form 

u — = 

for C2 = 4c'J, Sc^J = 3. 



(200) 

(201) 
(202) 



Painleve I as compatibility equation Inspired by the Saff-Totik approach, we construct 
the wavefunctions based on monic polynomials, (Pol is the polynomial part) 

n-l 

0n = n '^"(■^) = + Oiz""-'), (203) 



i=0 



and rewrite the equations for the Lax pair as 



-^0n = 0n+l + 



^n-2'"n-l 



0n+2 



2 T/i— — r " -n— IT'S— J- ' 2 

Notice that using the shift operator W, the system can also be written 

Introduce the scaling function through 

(pniz) = e2ft'0(C,z/). 



(204) 
(205) 
(206) 
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The action of Lax operators on t/j gives the representation 



- + h^x. = z + hdc = - + + h^^dr. 

2 ^ 2 ^ 



Therefore, the action of C, is given by the sum of equations at order /i^": 

- . . Ir^ .VV ^''^ 

2 2 
and the action of 9^ by their difference: 

,2 



Equivalently, we can write 

_ 1 
~ 2 

1 

Expanding the shift operator in h leads to 



(>V + lf +(rLi>V-^ + l)' 

(>v - 1)' - K-i w-^ - 1) 



4, 



Zi D 



and 



Substituting into the equations for dc^ gives the system of equations 

where primed variables are differentiated with respect to The equations can be written 
in matrix form as 



(207) 

(208) 

(209) 

(210) 
(211) 

(212) 
(213) 

(214) 



where 



A 



2(; + »j(2(,-M) 
9 



2C-M 
3 

u 

6 





2(C+") 



1 





The compatibility equations 

A - g' = [g, A] 

yield the Painleve equation derived in the previous section: 





and 



A 



[Q, A] 



+ 



6 



u 
6 

2(C+m)m 



" 3 
u 
' 6 



? 

3 ^ 



(215) 



(216) 



(217) 



(218) 



(219) 
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Thus 



= A - g' - [g, A] = ( 




2 
3 











) 



(220) 



V 6 9 

The only non-trivial element of the matrix gives 
u —Auii — 4 = 0, 



(221) 



i.e. the Painleve equation derived in the previous section. 

6.1.3. Conclusions The derivations presented above indicate that, in the vicinity of a 
(2, 3) cusp, the refined asymptotics for Laplacian growth are based on the behavior of 
the Baker- Akhiezer function for the Painleve 1 equation. This fact allows to properly 
define the evolution of the domain beyond the critical time, by identifying the support 
of the measure with the support of the zeros of this function. This is a work in progress 
which will be reported elsewhere. 

It is also interesting to note that the double-scaling limit required to derive the 
refined asymptote mirrors an earlier result, due in its original form to Stahl |150] . and 
related to orthogonal polynomials in |151j . It describes an approximation of the Cauchy 
transform of a planar domain via a special sequence of Fade approximants (in the spirit 
of section 5.3.2), which by exponentiation would translate into the double-scaling limit 
presented in this section. 

6.2. Non- commutative probability theory and 2D quantum models 

We conclude this review with a brief presentation of outstanding problems in two- 
dimensional quantum models, where the use of random matrix theory led to important 
results, and (perhaps most importantly) pointed out to the need for a probability theory 
for non-commutative random variables. In turn, such a theory is intimately related to 
the semi-normal operator approach presented in the previous chapter. 

6.2.1. Metal-insulator transition in two dimensions The details of the transition from 
conductive to insulating behavior for a system of interacting 2D electrons, in the presence 
of disorder, referred to as metal-insulator transition (or MIT) are not well understood, 
despite decades of research. Here we give a very sketchy description of this problem, in 
order to illustrate the mathematical essence of the model and of the difficulties, and we 
refer the reader to one of the several excellent monographs on the subject |152] . The fact 
that a system of electrons may "jam", i.e. behave like an insulator, because of either 
strong interactions (Mott transition) or strong disorder (Anderson transition) has been 
known for roughly half a century. However, creating a theoretical model which could 
incorporate both interactions and disorder in a proper fashion, was difficult to achieve. 
The foundation for our current formulation of this problem was laid by Wegner [H], 
and later improved by Efetov [T3]. A very clear exposition of this formulation can be 
found in the synopsis [16]. 
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In its simplest formulation, the model consists of a lattice in rf— dimensions (which 
may be taken to be Z*^), where to each vertex corresponds an ra— dimensional vector 
space of states (also called orbitals), and with hamiltonian 

H = Ho + Hd, Ho = ^ t^,Jx,n){y,n\, = ^ P\x,i){xJ\, 

n,{x,y) x,i,j 

where the state \n,x) depends on position x and orbital n, Hq refers to the interaction 
between adjacent vertices and implements the disorder component, via the random 
matrix which can be Hermitian, Orthogonal, Symplectic, etc. based on symmetries 
of the system. Efetov's idea was to use supersymmetry to incorporate interactions 
and disorder on the same footing; the method was later extended to implement 
the "Hermitization" of non-Hermitian random matrices with non-Gaussian weights, 
appearing in the same physical context |153j . For rotationally invariant measures, the 
authors showed that the distribution of eigenvalues can be either a disc or an annulus, 
and that there is a phase transition between the two, as a function of model parameters. 

The difficulties related to this formulation of the problem are due to the fact that 
the transition cannot be described within the established models of phase transitions. 
In all these models, the state of the system is obtained by minimization of a proper 
thermodynamic potential (for instance, free energy), or equivalently, finding the points 
of extrema of action in a path-integral approach (via a saddle-point condition). "Proper" 
phase transitions are characterized by potentials that are globally convex, so that the 
minimization problem is well-defined. However, the supersymmetric formulation of MIT 
does not lead to a true extremum, but rather a saddle-point, due to the non-compact, 
hyperbolic geometry structure of the effective theory {SU (1, 1) in the simplest case). The 
interested reader can find a detailed exposition of this phenomenon in |154] . for example. 
In the case where the system has a finite scale, it can be shown, following Efetov, that 
only the zero modes of the theory are important, which leads to an effective simplification 
in computing the multipoint correlation functions. However, the full model is still not 
solved for the 2D case, due to the difficulties pointed out above. 

In a nutshell, we may summarize the problem as non-tractable using the standard 
statistical physics formulation of phase transitions. In that sense, the situation is similar 
to another famous unsolved physical model, the disordered spin problem in the presence 
of magnetic field, where determination of the ground state is a task of exponential 
complexity (with respect to the size of the system). The phase transition where the 
system goes from an ordered state to a state with local order but no long-range order 
(a spin glass) is equally intractable as MIT, for the reasons explained. 

Interestingly enough, both problems may be approximately studied using a 
physicist's approach notorious for its lack of control: the replica-symmetry breaking 
(RSB) |155] . We mention it here mainly because of its statistical interpretation. 

Starting from the elementary observation logZ = lim„_^o(^'^ ~ is tempting 

to replace averages (over disorder) computed from the thermodynamic potential (log Z), 
with averages computed with (2'"'), because of the implicit assumption that repeated 
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products of the random variable Z will self-average (an implicit application of the 
Central Limit Theorem). By extension, one may assume that averages of products 
of operators, projected on special states, (O|0i02 • • • 0fc|O) (correlation functions), may 
also be computed by the same argument. 

In this (statistical inference) approach, the failure of standard descriptions of phase 
transitions is related to reducing correlation functions of products of operators, to their 
projections onto selected states. At the critical point, such projections do not have the 
expected convergence properties. It is therefore natural to ask whether one may use 
other (weaker) criteria to determine the critical point. In particular, is it possible to 
define statistical inference for the operators themselves, rather than special projections? 

The answer is affirmative, and such a theory was constructed almost in parallel 
with the MIT and spin glass models described above. 

6.2.2. Non- commutative probability theory and free random variables The basic 
elements in the probability theory for non- commutative operators |156j are the following: 
A, a non-commutative (operator) algebra over C, 1 G a functional (p : A ^ C, 0(1) = 
1, called expectation functional. 

Quantum mechanics offers specific examples: 

• Example 1: ^ = bounded operators over Hilbert space of states H, ^ G 7Y, ||^|| = 1, 
the ground state, and 

<f>{A) = m\o. 

• Example 2: A = von Neumann algebra over 7i, and functional = Tr. 

In order to develop inference methods within this theory, it is necessary to define 
the equivalent of independent variables in commutative probability. Such variables are 
called free, and satisfy the following property: Ai, A2, . . . , are free if 4>{Ai) = 0, and 

0(^11^12 ■ ■ ■ ^ifc) = 0, Ai- 7^ Aj+i- 

Using these tools, generalizations of standard results in large sample theory are 
possible. We mention a few: 

• The "Gaussian" distribution (limit distribution for Central Limit Theorem) in free 
probability theory is given by operators with eigenvalues obeying the semi-circle 
distribution (Wigner-Dyson) p(A) = ^ 0? — A^; 

• Similarly, the Poisson distribution has as free correspondent the operators with 
eigenvalues distribution according to the Marchenko-Pastur (elliptical law), p(A) = 
^(A-a)(6-A); 

• The free Cauchy distribution is the Cauchy distribution itself. 

Likewise, there is a notion of free Fisher entropy, Cramer-Rao bound, etc. 

As announced earlier, the relation between this theory, random matrices and 
operator theory for 2D spectral support is two-fold: on one hand, we have the 
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important result that random matrices, in the large size limit, become free non- 
commutative random variables. Thus, inference in free non-commutative probability 
may be approximated using ensembles of random matrices, which explains the success 
of this concept in the physics of disordered quantum systems. 

On the other hand, the limit distributions specified above (via Wigner-Dyson, 
Marchenko-Pastur laws, and their 2D counterparts), are described through spectral 
data. Taking this as a starting point, it is relevant to construct sequences of operators 
which approximate the spectrum, which points directly to the methods of section 5. 

As a last remark, an early attempt to employ non-commutative probability theory 
in MIT was reported in |157j . It is likely that the application of this generalized inference 
method will help elucidate open questions like the ones discussed in this section. 
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